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Kurzfassung 


Sheet Molding Compounds (SMC) sind diskontinuierlich faserverstarkte Ver- 
bundwerkstoffe, die aufgrund ihrer Fahigkeit, Verbundbauteile mit langen 
Fasern zu geringen Kosten zu realisieren, weit verbreitet sind. Sie ermóglichen 
Funktionsintegration, wie etwa den Einsatz von Rippen oder metallischen 
Einsátzen, und kónnen mit kontinuierlichen Kohlenstofffasern gemeinsam 
verarbeitet werden, um die Formbarkeit von SMC mit den überlegenen 
mechanischen Eigenschaften von kontinuierlichen Fasern zu kombinieren. 
Das Streben nach hochintegrierten und komplexeren SMC-Bauteilen er- 
fordert jedoch ein tiefes Verständnis der Verarbeitungsmechanismen und 
deren Einfluss auf die Leistungsfahigkeit eines Bauteils. Prozesssimulationen 
adressieren diesen Punkt, indem sie mógliche Fertigungsfehler und Prozess- 
parameter vorhersagen. Diese Ergebnisse kónnen nicht nur zur Prozessausle- 
gung und zur Reduzierung von Trial-and-Error-Phasen genutzt werden, son- 
dern auch für die anschließende Struktursimulation durch eine virtuelle 
Prozesskette. 


In dieser Arbeit wird die Prozesssimulation von SMC zunächst mit einem 
makroskopischen Referenzmodell auf Basis von Faserorientierungstensoren 
adressiert. Dies entspricht dem Stand der Forschung, aber die zugrun- 
deliegenden Annahmen von geraden Fasern, die viel kürzer als jedes ge- 
ometrische Merkmal sind, werden in anspruchsvollen SMC-Anwendungen 
oft verletzt. Dies führt zu der Hypothese, dass eine direkte Simulation einzel- 
ner Faserbündel erforderlich ist, um den SMC-Formfüllprozess komplexer 
Geometrien genau zu beschreiben. Basierend auf dieser Hypothese wird 
eine neuartige direkte Bündelsimulationsmethode (DBS) vorgeschlagen, die 
eine direkte Simulation auf Komponentenebene ermóglicht und dabei die 
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Beobachtung nutzt, dass Faserbiindel während des SMC Fließpressens oft 
in einer gebiindelten Konfiguration verbleiben. Das entwickelte DBS Mod- 
ell kann mit Patches kombiniert werden, um den Co-Molding-Prozess von 
SMC mit kontinuierlichen Faserverstárkungen zu simulieren. Daher wird 
ein Modell zur Beschreibung des Materialverhaltens von unidirektionalen 
Kohlenstofffaser-Patches einschließlich eines einfachen Schádigungsmodells 
zur Vorhersage von Defekten entwickelt. 


Die Parameter des makroskopischen Referenzmodells, des DBS-Modells und 
des Patch-Modells werden experimentell bestimmt. Dazu gehóren die ther- 
mischen Eigenschaften des SMCs, die temperaturabhängige und ratenab- 
hángige Viskosität der SMC-Paste, die Reibung an der Werkzeugwand sowie 
die Kompressibilitat des SMCs. Ebenso werden die temperaturabhángigen 
und ratenabhängigen mechanischen Eigenschaften der Patches bestimmt, 
die jedoch große Streuungen zwischen den Proben und Chargen aufweisen. 


Schließlich werden die Modelle auf mehrere Validierungsfälle angewandt, 
um die Anwendbarkeit auf Komponentenebene zu bewerten. Die Beispiele 
zeigen eine gegenüber dem makroskopische Referenzmodell verbesserte 
Vorhersage der Faserarchitektur, insbesondere der Faserorientierung in der 
Nähe von Werkzeugwänden sowie der Vorhersage von Bindenähten und 
Fließmarken. Zusätzlich bietet das DBS Modell die Option, Krümmungen der 
Bündel vorherzusagen und den Faservolumenanteil zu berechnen, welche 
durch Mikro-Computertomographie, thermisch gravimetrische Analysen 
und Durchleuchtungsbilder validiert werden. 
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Abstract 


Sheet Molding Compounds (SMC) are discontinuous fiber reinforced compos- 
ites that are widely applied due to their ability to realize composite parts with 
long fibers at low cost. They enable function integration, e.g. ribs or metallic 
inserts and can be co-molded with continuous carbon fibers to combine the 
formability of SMC with the superior mechanical properties of continuous 
fibers. However, the pursuit for highly integrated and more complex SMC 
components requires a profound understanding of the processing mech- 
anisms and their influence on the performance of a component. Process 
simulations address this point by predicting possible manufacturing defects 
and process parameters. These results can not only be used to configure the 
process and reduce trial and error phases, but also for subsequent structural 
simulation by virtue of a virtual process chain. 


In this work, the process simulation of SMC is addressed initially with a 
macroscopic reference model based on fiber orientation tensors. This is in 
line with the state of research, but the underlying assumptions of straight 
fibers that are much shorter than any geometric feature are often violated 
in advanced SMC applications. This leads to the hypothesis that a direct 
simulation of individual fiber bundles is required to accurately describe the 
SMC mold filling process of complex geometries. Based on this hypothesis, a 
novel Direct Bundle Simulation (DBS) method is proposed to enable a direct 
simulation at component scale utilizing the observation that fiber bundles 
often remain in a bundled configuration during SMC compression molding. 
The developed DBS model can be combined with patches to simulate the 
co-molding process of SMC with continuous fiber reinforcements. Hence, 
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a model is developed to describe the material behavior of unidirectional 
carbon fiber patches including a damage model to predict defects. 


The parameters of the macroscopic reference model, DBS model and patch 
model are determined experimentally. This includes thermal properties of 
the SMC, temperature dependent and rate dependent viscosity of the paste, 
friction at the mold surface as well as the compressibility of the SMC. Likewise, 
the temperature dependent and rate dependent mechanical properties of the 
patches are determined, but they show large scatter between samples and 
batches. 


Finally, the models are applied to several validation cases for evaluating the 
applicability at component scale. The examples show improved prediction of 
fiber architecture compared to the macroscopic reference model, especially 
the fiber orientation near mold walls and the prediction of knit lines or flow 
marks. In addition, the DBS model provides the option to predict bundle 
curvatures and calculate fiber volume fractions, which are validated by micro- 
computed tomography, thermal gravimetric analysis and fluoroscopy. 
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1 Introduction and motivation 


Lightweight design is a holistic development approach, which ideally triggers 
a positive self-reinforcing effect: If mass is saved, inertial loads on a structure 
and loads from powering a machine are reduced. This scales down require- 
ments on support structures, which then can be designed lighter themselves 
to reduce mass even further. Ultimately, this reduction in mass reduces 
primary energy consumption for all mass related loads in a machine and 
decreases impact energy in a catastrophic failure. 


One option to reduce mass is the application of composite materials. Such 
materials combine properties of multiple materials to obtain a tailored prop- 
erty profile suited to the specific application. Fiber reinforced plastics (FRP) 
are an important group of composites comprising strong load-carrying fibers 
(typically made from glass or carbon) and a matrix material (thermoset or 
thermoplastic polymer), that bonds fibers together, transfers loads to fibers 
and protects fibers from environmental influences [1]. If such reinforcement 
fibers span the entire length of a part, the reinforcement type is classified 
as continuous FRP (CoFRP). Continuous fiber reinforcements offer high po- 
tential for mechanical performance, especially regarding strength, due to 
large fiber volume contents and high alignment of fibers. Consequently, they 
are a popular choice in the aerospace industry, performance cars and sports 
products requiring high stiffness. However, the necessary raw materials of 
CoFRPs are rather expensive, the degree of process automation is limited, 
and careful consideration is needed to make economically viable material 
choices [2]. Additionally, the geometric shapes are restricted by the underly- 
ing architecture of fabrics - sometimes it is simply not possible to place fibers 
such that they span an entire part. Chopped fibers, on the other hand, offer 


1 Introduction and motivation 


much more design freedom as they can be molded together with the matrix 
material in a flow process such as injection molding or compression molding. 
Such a material is classified as discontinuous FRP (DiCoFRP). DiCoFRPs can 
be molded to complex geometries at high rates and thus are very cost efficient. 
However, their mechanical properties are inferior to CoFRPs. Discontinuous 
fiber reinforced plastics with local continuous fiber reinforcements (CoDi- 
CoFRP) are a new class of composites aiming to combine the merits of both 
constituents [3]. They use continuous high performance fibers in critical 
loading areas and a discontinuous compound to realize function integration 
with ribs, beads or other complex geometric features. 


The mechanical properties of FRPs depend significantly on the manufactur- 
ing process, which affects key properties such as fiber orientation and fiber 
volume fraction. Prediction of these properties by numerical simulations 
can help to account for these effects in an early product development phase 
and to save expensive trial and error studies. Process simulations aim to 
predict manufacturing defects and material properties after processing an 
FRP Ideally, these results are fed to subsequent structural simulations in the 
framework of a continuous CAE chain for more accurate results of structural 
simulations [4-6]. 


Therefore, this work addresses the simulation of a Sheet Molding Compound 
(SMC) CoDiCoFRP manufacturing process. Most current models for SMC 
compression molding simulations assume a macroscopic behavior without 
accounting for the underlying fiber bundle architecture of SMC at mesoscale, 
which can lead to inaccurate predictions of fiber orientations after processing 
in complex geometries. The objective of this work is an accurate prediction of 
the manufactured fiber bundle architecture and the possibility to incorporate 
local unidirectional reinforcement patches in the simulation. The developed 
methods should result in a better understanding of defects during the early 
design phase of SMC components and thus pave the way to more structural 
applications of CoDiCoFRP. These applications have the potential to reduce 
primary energy consumption in a cost-effective manner. 


2 Fundamentals and state of 
research 


2.1 Modeling scales 


Fiber reinforced composites can be analyzed and modeled at different scales 
ranging from the molecular scale at which chemical processes occur (curing, 
fiber sizing) up to the component scale at which a part is loaded. This work 
terms the smallest scale microscale, which refers to the scale of individual 
fiber filaments with a diameter of approximately 15 um. The term mesoscale 
refers to fiber bundles or fiber tows that consist of hundreds of individual 
fiber filaments. The dimensions of fiber bundles cross sections are in the 
order of 0.1 mm, while the length is 25 mm to 50 mm. Modeling approaches 
at mesoscale treat the entire bundle as a single instance with no resolution 
of the underlying micro-structure. The macroscale refers to a homogeneous 
macroscopic scale, at which the micro- and meso-structure of the composite 
is completely blurred and described by effective macroscopic properties. The 
scales in this work are summarized in Figure 2.1. 


Micro Meso Macro 


Figure 2.1: Schematic illustration of the modeling scales in fiber reinforced composites. Fibers 
are at the micro scale, bundles at the mesoscale and the homogeneous material with blurred 
micro- and meso-structure is a representation at the macroscale. 


2 Fundamentals and state of research 


2.2 Fiber suspensions 


This section starts with an introduction to general properties of fiber suspen- 
sions and subsequently states deformation measures and balance equations. 
After introducing stochastic descriptions of the fiber orientation state at 
macroscale, the section summarizes approaches to model rheology of fiber 
suspensions at macroscale. Finally, direct numerical models are reviewed 
and classified. 


2.2.1 General properties and classification 


A fiber suspension is a material consisting of elongated particles suspended 
in a viscous matrix material. Such particles are characterized by their aspect 
ratio 
fp -— (2.1) 
P q 


with the (equivalent) diameter d and length L. The aspect ratio describes the 
dimensionless slenderness of a particle. Another dimensionless shape factor 
is defined as 


mcd 
£5 Q.2) 
rtl 


and simplifies the notation of fiber orientation evolution equations later in 
this section. 


A suspension of reference volume V; that contains Np particles can be de- 
scribed either by the number density 


=e (2.3) 
Np = — : 
P y 
or the fiber volume fraction 
nd2L 
f= Vpmp = 4 "P (2.4) 
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where Vj described the volume occupied by one cylindrical fiber. 


Fiber suspensions are commonly divided in different concentration regimes 
[7] and different theories apply in each of these regimes: 


Dilute regime Particles are spaced at a large distance compared to their 
dimensions. Thus, the particles can be treated individually, regardless 
of the existence of other particles. The effect on the surrounding fluid 
flow is typically neglected. The dilute concentration regime is formally 
bound by 


1 

Np K I (2.5) 

Semi-dilute regime Particles interact with each other through flow field per- 
turbations. They can be distributed isotropic and are well below the 
nematic transition [7]. The semi-dilute regime is formally bound by 


1 1 
B « np K PTA (2.6) 


Concentrated regime The average distance between particles is of the order 
of their diameter. They interact via hydrodynamic effects and contacts 
[8]. The concentrated regime is formally bound by 


1 
22 < Np. (2.7) 
Typical industrial composites belong to the concentrated regime, as 
illustrated in Figure 2.2. 


Nematic regime (liquid crystal) Particles are suspended freely but remain 
in a constrained alignment, as this is the only possibility to achieve a 
corresponding high fiber volume fraction. 
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Figure 2.2: Concentration regimes: The black lines separate the dilute regime (white) from the 
semi-dilute regime (light gray) and the concentrated regime (gray). 


The bending stiffness of a flexible cylindrical fiber in a suspension may be 
described by the dimensionless ratio 


En 


S* = —— 
Anyro 


(2.8) 
with the Young's modulus of a fiber E, scalar shear rate y and matrix shear 
viscosity n [9]. Fibers bend with an increased aspect ratio, higher shear rates 


and reduced bending stiffness [10], which is summarized by a small value of 
S*. 


2.2.2 Rate of deformation measures 


The motion of a fiber suspension is characterized by its differentiable velocity 
field v : D — R? on a domain D c R?. The position within that domain is 
denoted by x € D and differential operators work with respect to æ, if not 
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mentioned otherwise. Relevant objective measures for the deformation are 
the symmetric strain rate tensor 


1 
D= 5 (gradu + gradu" | (2.9) 
and the vorticity tensor 
1 T 
W = A (gradv — grado » (2.10) 


where grad denotes the gradient operator. The spherical part of the strain 
rate tensor i 
D'-;u(D)I-P:D (2.11) 


describes the volumetric rate of change and vanishes in case of incompress- 
ible fluids. The tensor P; = i & J denotes the spherical projection tensor of 
fourth order. The deviatoric part of the strain rate tensor is denoted as 


D'D-D'z:D (2.12) 


and describes the rate at which the fluids shape changes. The tensor P2 = 
I5 — P; denotes the deviatoric projection tensor of fourth order and Ii? is the 
symmetric part of the fourth order identity tensor. A scalar measure for this 
deformation rate is the shear rate y = V2D' : D'. 


2.2.3 Balance equations 


The solution fields of the deformation of fiber suspensions are described by 
several balance equations. These balance equations describe conservation 
properties observed in nature, e.g. conservation of mass, momentum and en- 
ergy. A general form of such a balance equation for an arbitrary differentiable 
scalar field (e) : D — R in an Eulerian frame is 


ö(e) , z 
ES t div((e)v) = s (2.13) 
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with the spatial time derivative os) at a fixed point in space, the divergence 


operator div(e) and source term s. The source term s summarizes the pro- 
duction of the property, non-convective flux of the property and externally 
applied sources. The spatial time derivative is related to the material deriva- 
tive via 
. _ Oe) 
(e) = ET +v: grad(e) (2.14) 


and describes the change of a property at a deformed material point. 


The first important field of interest is the mass density p : D — R. As mass can 
be neither produced (at least not in Newtonian mechanics) nor transported 
in other ways than by convection, the source term vanishes and leads to 


ô 
E 4 div(pv) = 0 (2.15) 


for the balance of mass in an Eulerian frame. Application of Equation (2.14) 
gives the Lagrangian form 
p 4 pdiv(v) 2 0 (2.16) 


with the material derivative of the mass density field p. 


Another balanced property is the momentum density field pv and using the 
tensor form of (2.13) results in the balance equation for momentum 


“ee + div((pv) 9 v) = div(o) + pb (2.17) 


with the non-convective momentum flux density o also known as Cauchy 
stress tensor and a body force density field pb. Using the mass balance (2.15) 
and (2.14), this equation may be recast to the Lagrangian form 


p? = div(o) + pb. (2.18) 
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Finally, the temperature field T : D — R is of interest and described by the bal- 
ance of internal energy pcp T, where cp is the specific heat capacity. Without 
heat sources, this is 


O(pcp T) 


* div((pcy T) -v) = -grad(d) - o : D (2.19) 


with the non-convective heat flux d, and the dissipated strain energy o : D. 
For a constant specific heat capacity, the Lagrangian form of Equation (2.19) 
is 

pcp T 2 - div(d) e : D. (2.20) 


2.2.4 Stochastic description of fiber orientation states 


The orientation of a single straight fiber, fiber segment or bundle segment 
is denoted as a unit vector p € S, where S = (pe R° | 1 = |||]. The change 
of fiber orientation with time p is of upmost interest for process simulations 
of fiber suspensions. The first model to describe the orientation evolution 
of a single rigid spheroid suspended in an infinitely sized Newtonian fluid 
domain without buoyancy and inertia was described by Jeffery in 1922 as 


py =W-p+é(D-p-(D:pep)-p| (2.21) 


with the symmetric strain rate tensor D and the skew-symmetric vorticity 
tensor W of the suspending fluid [11]. For pure shear in a 2D domain, this 
equation reduces to a scalar evolution equation for the fiber orientation angle 
0 


Be r (1+ £cos20). (2.22) 


Strictly speaking, Jeffery’s model is only valid for a single ellipsoid in an 
infinite domain. However, it can be shown, that the equation is also valid for 
cylindrical fibers, if an equivalent aspect ratio re is used [12]. Cox [13] derived 
such an equivalent aspect ratio based on slender-body theory. Zhang et al. 
[14] utilized Finite Element Analysis to propose a fitted cubic relation between 
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re and rp, which is more accurate for small aspect ratios. Jeffery’s model yields 
reasonable results in the dilute regime, where fibers are dispersed so sparsely 
that they do not affect each others motion. 


Modeling the orientation evolution in semi-dilute fiber suspensions requires 
a formal treatment of an orientation state consisting of several fiber directions. 
This can be achieved either by a stochastic description, as described in the 
following, or modeling of individual fiber instances, as described in the later 
Section 2.2.6. 


The fiber orientation density distribution function V : S — P gives the proba- 
bility to find a fiber in direction p with P = (V € R| 0 < V < 1}. It fulfills the 
normalization condition 

[voa =1 (2.23) 


and the continuity condition 
Y=-grads(P(p)D). (2.24) 


Here, dp is the surface element on the unit sphere S that ensures an invariant 
integration and grad, is a gradient on the unit sphere S. Further, V(p) is a 
symmetric function, i.e. Y(p) = Y(-p), as the end point of a fiber cannot be 
distinguished from its starting point. 


Folgar and Tucker modified Equation (2.21) to 


V = — grads (V (p)pj) + Giy grads (V (p)) (2.25) 


by employing the fiber orientation density distribution function and adding 
a term C;Y for diffusion [15, 16]. The parameter Cj is termed the interaction 
coefficient and for C; = 0, Equation (2.25) is identical to Equation (2.21). The 
diffusion term is inspired by Brownian particle motion and models the ten- 
dency of fibers to distribute to an isotropic state due to interactions between 
individual fibers. The empirical parameter Cj describes the rate at which 
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fibers move towards a more isotropic state and depends on fiber aspect ratio, 
fiber volume fraction and other properties of the fiber suspension. [15]. 


The fiber orientation density distribution function V (p) may have arbitrary 
complexity. Thus, it is inconvenient to store in numerical simulations and 
Advani and Tucker [17] suggested to store moments instead. These moments 
are called fiber orientation tensors and commonly used moments are the 
second order fiber orientation tensor 


A= | v@pepdp (2.26) 
and fourth order fiber orientation tensor 
A= | vppepepe pap. (2.27) 
Fiber orientation tensors are fully symmetric and the contractions 
A-I=1 A:I=A (2.28) 


apply. The set of second order orientation tensors fulfilling such properties is 
termed A. Exemplary fiber orientation states, corresponding second order 
fiber orientation tensors and plots of the fiber orientation density distribution 
function are given in Figure 2.3. 


Jeffery’s model (2.21) and the Folgar-Tucker model (2.25) can be expressed in 
terms of fiber orientation tensors as 


A=W-A-A-W+é(D-A+A-D-2A:D) (2.29) 
and 
A=W-A-A-W+é(D-A+A-D-2A:D)+2CGy(I-3A), | (2.30) 


respectively. However, the fourth order tensor A is required to solve these 
equations. One could obtain it from analogous equations for À, but that 
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1 1 
2 : 0 3 i 0 
A-|0 5 0 A-|0 3 : 
0 0 0 0 0 1 
(d) Unidirectional tensor (e) Planar tensor (f) Isotropic tensor 


Figure 2.3: Orientation examples. The distributions (a), (b) and (c) are generated by the procedure 
described in Section 5.1.2. 


would require a sixth order tensor and so forth. Hence, A is typically obtained 
from a closure approximation. The simplest analytical closure is a quadratic 
closure 

A9 - Aa A, (2.31) 


which is correct only for unidirectional fiber orientation. The linear closure 


AL- 


(Ae1-19A«A I-IDA4 (ADD + (IDA)? 


Sir 


(2.32) 


E IeI-IUI-*UUDI'* 
35 


is exact for a random isotropic fiber orientation state. Hybrid closures try to 
interpolate these states [17, 18]. More advanced closures approximate the 
fourth order tensor by fitting eigenvalues such as orthotropic fitted closures 
[19, 20] and the invariant-based optimal fitting closure (IBOF) [21]. This work 
proceeds to use the IBOF closure in macroscopic reference solutions, as it is 
a good compromise between accuracy and computational cost. 


The Folgar-Tucker model was extended by several authors and a list (non- 
comprehensive) is given in Table 2.1. In non-dilute suspensions, Equations 
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(2.29) and (2.30) overestimate the reorientation rate of fibers in comparison to 
experimental evidence. Thus, variants with reduced strain rate (RSC) or a re- 
tarded principal rate (RPR) have been developed and add up to two additional 
parameters. Additionally, the isotropic diffusion term has been replaced by 
an anisotropic term [22], which introduced four additional parameters. A 
Maier-Saupe potential was suggested to counteract the diffusive process 
considering nematic conditions in suspensions with high fiber volume frac- 
tions [23]. There are also some suggestions to account for confinement of 
fibers in regions, where the distance between mold walls is smaller than the 
length of suspended fibers [24, 25]. Recent developments tend to reduce the 
number of parameters for easier parameter identification [26, 27]. 


Most of these models are developed primarily for injection molding simu- 
lations with a three-dimensional short fiber architecture in mind and the 
results can be tuned by the choice of parameters. To obtain an unbiased ref- 
erence solution, the macroscopic models later in this work use Jeffery's basic 
equation without retarded principal rates and without diffusion parameters. 


Table 2.1: Overview on stochastic macroscopic fiber orientation models and number of required 
parameters for the diffusion model and the strain rate reduction. Implementations and examples 
may be found at https: //github.com/nilsmeyerkit/fiberoripy 


Model Year Diffusion Strain rate Reference 
Jeffery 1922 0 0 11 
Folgar-Tucker 1984 1 0 15 
RSC 2008 1 1 28 
ARD-RSC 2009 5 1 22 
FTMS 2010 2 0 23 
iARD-RPR 2016 2 2 29 
pARD-RPR 2017 2 1 26 
MRD 2018 3 0 27 


Stochastic fiber orientation descriptions are valid for fibers that are much 
shorter than any dimension of the flow domain, as they do not account for 
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confinement by mold walls and other non-local effects. A sufficiently high 
number of fibers is needed at each point, at which an evolution equation 
is evaluated, to represent a proper sample of the fiber orientation density 
distribution function. Additionally, these models typically assume a constant 
fiber volume fraction and constant fiber aspect ratio. 


2.2.5 Rheology of fiber suspensions 


The total macroscopic stress in an incompressible fiber suspension may be 
decomposed into 
o=—plt+oytomt op: (2.33) 


with a contribution from the volumetric pressure p, from the neat matrix 
fluid eun a deviatoric extra stress due to fiber-matrix interactions TiM (long- 
range interactions) and a deviatoric extra stress due to fiber-fiber interactions 
Ou (short-range interactions, such as mechanical contacts and lubrication) 
[30,31]. The illustration in Figure 2.4 visualizes the sources of long-range 
hydrodynamic interactions due to disturbances of the flow field and short- 
range interactions at contact points of fiber bundles. 


Two-way coupled 
hydrodynamic interaction 
(long-range or fiber-matrix) 
Contact, friction, lubrication 

(short-range or fiber-fiber) 


Figure 2.4: Illustration of long-range hydrodynamic interactions of bundles with the matrix 
velocity field and short-range fiber-fiber interactions at contact points between bundles. The 
hydrodynamic interaction is two-way coupled, i.e. the motion of a bundle depends on the 
fluid velocity in its environment and the fluid velocity is affected by bundle motion. Adopted 
from [32]. 
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2.2.5.1 Matrix contribution 


The constitutive equation for the neat matrix fluid 


om = 6tr(D)I +2nD' = (37P1 + 2nP2): D (2.34) 
Viso 

computes oy from a scalar shear viscosity 7 and a scalar bulk viscosity ¢, 
which may be combined to a fourth order isotropic viscosity tensor Viso. The 
term tr(D)I = P, D vanishes in incompressible fluids and most compressible 
engineering cases assume ¢ = 0, as ¢ is only relevant if the reciprocal of the 
bulk strain rate was at the order of magnitude of molecular movement [33]. 
Therefore, Equation (2.34) is reduced to 


oy =2qD' (2.35) 


for the following sections. Subjecting the isotropic matrix to uniaxial, incom- 
pressible elongation (Dj, = D, = -D/,./2) results in the first and second 
normal stress differences 


Nm = OMxx 7 My = 30D. (2.36) 
NM2 = 04x zz = 38ND xx (2.37) 


of the matrix. The factor 37 may be termed elongational viscosity of the 
isotropic matrix 7) y, and the ratio between this elongational viscosity and 
the shear viscosity 7) m,xx/n is known as Trouton ratio, which is equal to 3 for 
isotropic incompressible Newtonian media [33, 34]. 


2.2.5.2 Long-range contribution due to hydrodynamic interaction 
of fibers with the matrix 


Typical fibers in composite applications have a high Young's modulus and 
thus are considered inextensible compared to the suspending matrix material. 
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Hence, fibers represent a constraint on the fluid field, disturbing the velocity 
field at fiber scale and leading to additional shear deformation compared 
to the macroscopically applied deformation rate. This extra shearing of the 
viscous matrix translates to extra stress at macroscale that depends on the 
fiber volume fraction, fiber aspect ratio and fiber orientation. This extra stress 
results purely from the long-range velocity disturbance even in absence of 
any mechanical contact of fibers. 


The general objective form of anisotropic stress in an incompressible fluid 
derived by Ericksen [35] and can be expressed as 


orm = —Pol +2n0oD+mAtnoA:D+2n3(A-D+D-A) (2.38) 


using orientation tensors [36]. The first two parameters Po and no describe 
the isotropic extra stress by fibers. The three parameters 7),72,73 are the 
internal fiber stress, elongational extra viscosity in fiber direction and dif- 
ference between shear viscosity along fibers and transverse to the fibers, 
respectively [37]. The term 7; is typically neglected under the assumption 
that there is no internal strain-rate-independent fiber stress at scales, where 
Brownian motion is insignificant. 


The extra stress from fiber-matrix interactions should not depend on the 
trace of D and should be trace-free itself, such that the term p describes the 
complete hydrostatic pressure [38]. Thus, Equation (2.38) is reformulated to 
the deviatoric stress 


2noP2+ n2 


ke Tei 
3 


2 
oly = +2ns[A I+I A-s1e4] ant. 


(2.39) 
This equation is obtained by replacing D with the deviatoric strain rate D’, 
multiplying Equation (2.38) with P» from the left and using the identity T : A = 
A. Forno = 0 and 73 = 0, i.e. slender fibers with negligible thickness in a dilute 
suspension, Equation (2.39) becomes equivalent to the equation derived by 
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Batchelor [39]. For dilute suspended elongated spheroids, Batchelor obtained 
the expression 
2f 


Ban (2.40) 
3 [imer - 1.5] 


n2 


using slender body theory [40-42]. Dinh and Amstrong suggest the variant 


2fr2 
I (2.41) 


3lnt/zif) 


for semi-dilute suspensions assuming aligned fibers [42, 43]. Shaqfeh and 


12, DA = 


Fredrickson [44] rigorously derived the expression 


afr 
n (2.42) 
3[In(1/ f) + InIn(1/ f) + C] 


12, SE = 


for slender fibers using multiple scattering, accounting for velocity distur- 
bances at multiple length scales. The factor C ranges from 0.1585 for aligned 
fibers to -0.6634 for isotropic fiber orientation [30, 44]. These authors also 
imply 79 = 0, since this contribution is small compared to 7». 


Tucker [41] recasts Equation (2.38) together with the isotropic matrix compo- 
nent to 


cw + aiu =2Ñ |D + NA: D N (AOD + DOA)| (2.43) 


with 7j containing the entire isotropic contribution from the matrix and sus- 
pended particles. The dimensionless particle number N, describes the im- 
portance of extra stress due to particle stretching and is usually dominant 
over the shear number Ng, which describes increased shear resistance due to 
the thickness of fibers [41]. Strictly enforcing deviatoric stresses and strain 
rate tensors leads to 

Ou + Op, = Va: D' (2.44) 
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with the fourth order anisotropic viscosity tensor 


Va = 27 | P2 + Np + Ng (2.45) 


À-l1faA 
3 


AOI +T A-=1eA| 


following Bertóti [38]. This version is equivalent to a combination of Equa- 
tions (2.39) and (2.35) for a proper choice of the parameters, but has a dif- 
ferent interpretation. Such parameters can be obtained e.g. from fitting 
micro-mechanical results [38,45] or by expressing them in terms of Equations 
(2.40) to (2.42). Equation (2.43) is also commonly used directly to model 
anisotropic viscosity with N; = 0 for slender fibers [23, 46]. 


Other authors do not start from the general form given in Equation (2.38), but 
with the transversely isotropic case of an unidirectional fiber reinforcement 
in analogy to linear elasticity [47, 48]. For an incompressible material, this is 


1 1 H 


B xx Nxx 2N xx 20x 000 O FM,xx 
Di, | ihe + m T m 0 0 0| [oru 
Da) -= ES i az) O 0 9 ||gmuz| 2.46) 
D yz T 0 0 OFM, yz 
Dix T 0 OFM,zx 
Diy symm, "m C EM,xy 


in Voigt notation and for fiber alignment in x-direction. Pipes et al. [48, 49] 
determined the corresponding parameters from micro-mechanical analysis 


to 
nf av f 2 
a= |— |r (2.47) 
2 \l-a/f 
for the elongational viscosity in fiber direction, 
n|[2-avf 
= | ———— (2.48) 
dE: G -a 2 
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for the axial shear viscosity and 


nyz=n(l-@f) (2.49) 


for the transverse shear viscosity. A geometrical factor a describes the ideal- 
ized arrangement of fibers (a? = V12/r for hexagonal packing and a? = 4/7 
for square packing). Similar to the volume averaging suggested by Advani and 
Tucker [17] for elastic stiffness, the result is then used in a volume-averaged 
equation that depends on fiber orientation tensors [50]. 


Asummary of the parameterizations for anisotropic viscosity is given in Table 
2.2. It is noteworthy that all variants rely on three parameters compared to 
the five parameters of transversely isotropic materials in elasticity due to 
the assumption of incompressibility. In principal, these formulations are 
equivalent to each other. The general form with 7,72,73 is the preferred 
parameterization in subsequent chapters, because it is based directly on the 
matrix viscosity, which can be measured in a simple rheometer. 


Table 2.2: Parameters for macroscopic viscosity including matrix and long-range contributions. 


Parameterization Equations Parameters References 
General form (2.35) and (2.39) N, N2, N3 [35, 39, 42, 44] 
Recasted general form (2.43) ñ, Np, Ns (23,41, 45, 46] 
UD + volume averaging (2.46) Nxx Nyz» Nxy [47, 48, 50,51] 


Figure 2.5 compares predicted first normal stress differences for a suspension 
of aligned fibers. Batchelor’s model and the Dinh-Amstrong model, which are 
designed for dilute and semi-dilute suspensions, underestimate the increase 
of elongational viscosity for higher volume fractions. Pipes’ model predicts a 
singularity as soon as the fiber volume fraction approaches the packing limit 
(square packing is assumed here). The Shaqfeh-Fredrickson model resembles 
the properties of Batchelor’s model for small volume fractions and those of 
Pipes’ model for larger volume fractions, including a singularity at the packing 


19 


2 Fundamentals and state of research 


4 10° 5 
z J 
ES j 
[5] - 
[=l 
z al 
& 1024 
=| = 
2 j 
3 | © Trouton ratio 
Dn " 
ir yita ——— Pipes 
= E Batchelor 
E F —— Dinh-Amstrong 
Es 4 Shaqfeh-Frederickson 
FH 
10° | | | | j j 


j j 
0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.45 
Fiber volume fraction 


Figure 2.5: Comparison of the first normal stress difference for a suspension of unidirectional 
aligned, square packed fibers with aspect ratio rp = 25. 


limit. It seems applicable to a relatively large range of fiber volume fractions 
and is therefore used in the reference model in this work. All models scale 
with n: which results in a significant increase of the elongational viscosity 
for high aspect ratios of the fibers. 


2.2.5.3 Short-range contribution due to fiber-fiber interactions 


Fibers in concentrated suspensions have multiple contact points at which 
fibers interact through mechanical friction or lubricated contact with a small 
matrix layer undergoing high shear. Such direct contacts between individual 
fibers are summarized as short-range interactions or fiber-fiber interactions. 
The macroscopic stress contribution can be formally expressed as 


1 M 
om=-7 Y, hap ® fap (2.50) 
T (a, p) 


with a gap vector hag and an interaction force fgg for all N. interaction pairs 
(a, B) of a reference volume V; [31,52]. The gap vector hap describes the 
shortest path between two straight fiber segment surfaces and is directed 


20 


2.2 Fiber suspensions 


from partner a to partner f. The scalar distance between two surfaces is de- 
noted g and may become negative, if fibers overlap. The interaction force fap 
comprises an elastic contact force in normal direction f and a tangential 
friction force fag as well as lubrication forces in normal direction er and 
tangential direction £s The focus of this section is the effect of averaged 


short-range interactions at a macroscopic scale. 


The average number of contact points per fiber is given as 


n Afri 01 (2.51) 

with orientation functions 
®| = Í INEST sin| Z (p,p rovo )dp'dp (2.52) 
$,- i [ cos(Z(p,p p’) V(p)V(p^dp'dp (2.53) 


as defined by Toll and Mánson [53-55]. Exemplary values of the orientation 
functions are given in Table 2.3. For a planar network of slender elastic fibers 
(rp > 1), the number of contact points per unit volume is consequently 


Ne 16f? 
— 877 2.54 
W^ we (2.54) 
The average normal force at contact points is then given as 
ME = S Edo T: (2.55) 


where E is the Young’s modulus of fibers and the averaging is indicated by a 
bar [53,54]. 


Servais et al. [56,57] build on these results for their investigation of short- 
range fiber-fiber interactions for dispersed fibers and bundled configurations. 
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Table 2.3: Exemplary values of the orientation functions [53]. 


91 ®2 
Unidirectional 0 1 
Planar isotropic 2/z 2/1 
Isotropic 71/4 1/2 


They propose a power-law relation for the tangential hydrodynamic lubrica- 
tion 


tl n-1 
fil, = ks | Avap| Avap (2.56) 


with power-law exponent n, hydrodynamic friction coefficient ks and the 
relative tangential sliding velocity Avgg. For the tangential friction, they 
suggest a Coulomb friction model 


sig =o sai] eo sn 


with a Coulomb friction coefficient u. The friction acts opposed to the di- 
rection of the sliding velocity and the operator [e] = (e)/|le|| is used as a 
convenient notation to compute the normal direction. 


Servais et al. [56] showed that Coulombic friction results in a yield stress 
of the macroscopic suspension with dispersed fibers and thus the result 
renders a Herschel-Bulkley fluid. For fiber bundles, they proposed a similar 
model using a Carreau relationship for the shear rate dependency. However, 
hydrodynamic lubrication dominates in case of fiber bundles, as the sheared 
surface area is large compared to the effects of the friction at the contact 
point [57]. 


A more general formulation of the short range stress contribution was ob- 
tained by Djalili-Moghaddam and Toll [30] for linear Newtonian lubrication 


forces 
3 


tl 4hy 
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with a constant interaction range 2hp, in which short-range interactions 
occur, and a parameter kp that accounts for all geometric details of the 
lubrication flow. Their final result 


4kp 2, 2m. 
OFF = 3,2 Pas :D (2.59) 


is independent of the interaction range [30]. The key term is the fourth order 
interaction tensor 


F [fle xp | Y(p)¥ (p')p® p® ps pdpdp' (2.60) 


that evaluates the interaction between fibers and vanishes for perfectly 
aligned fibers. Djalili- Moghaddam and Toll [30] compute the interaction 
tensor from discrete fiber directions (compare Section 2.2.6), but Férec et 
al. [58] suggest to approximately relate the interaction tensor to the known 
fiber orientation tensors. They denote the second order interaction tensor as 


B= Í Í |px p | YDY ppe papap (2.61) 
SJS 


and obtain the result " 
Bs 7 (A-A: A). (2.62) 


by replacing the Onsager potential | pxp | with the Maier-Saupe potential 
\|p x »'| ? and an appropriate prefactor [58]. They initially propose a quadratic 
closure 


o- lpeaB (2.63) 
$1 


and a linear closure [58]. The closures were later improved [59], but the effect 
of short-range interactions has been investigated in more detail with direct 
numerical models, as precise knowledge of the fiber architecture is necessary 
to properly account for this effect. 
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This work will use equations (2.59), (2.62) and (2.63) later in Chapter 4 for 
the macroscopic reference model. This set of equations has just a single ad- 
justable parameter kp and provides a closed formulation of the macroscopic 
stress contribution from fiber-fiber lubrication, which is rarely addressed in 
literature. 


2.2.6 Direct numerical models 


Direct numerical models solve the motion of individual fibers or representa- 
tive sets of fibers in a suspension instead of a stochastic moment. An equation 
of motion is solved for fiber segments, bead chains or other discretization 
forms of a fiber instead of an evolution equation for a fiber orientation tensor. 


2.2.6.1 Kinematic models 


The first class of direct simulation models transports fiber nodes according to 
the bulk velocity of the suspension. These models do not include the momen- 
tum balance or any forces and are therefore termed kinematic models here. 
One such model is commercialized as 3D TIMON Direct Fiber Simulation 
(DFS) feature by Toray Engineering Co., Ltd., Japan. A fiber is represented by 
a chain of rods connected with nodes as shown in Figure 2.6. These nodes 
are transported with a pre-computed velocity field vo, (x, t) of a macroscopic 
mold filling simulation in a forward Euler scheme as 


8r ca gs Ce (2.64) 
where ke N denotes a node index and n € N denotes a time step index [60]. A 
regularization is necessary to constrain fiber elongation and is achieved with 
a correction scheme ae > m based purely on the geometric configu- 
ration of the fiber [60]. This kinematic model has been used to analyze fiber 


bending in Long Fiber Thermoplastic (LFT) compression molding [61] and 
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shows good agreement of fiber orientation in compression molding of carbon 
prepreg platelets compared to Micro Computer Tomography (UCT) [62]. 
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Figure 2.6: The kinematic direct fiber simulation model computes a new temporary node posi- 


tion Br from a prescribed velocity field and the current node position xy by forward Euler 


integration. Afterwards, a regularization is applied to ensure constant fiber length. 


This model does not account for short-range interactions between fibers 
nor for long-range interactions by disturbance of the velocity field. This 
makes the model very efficient, but has the disadvantage that multiple fibers 
can occupy the same space and therefore leads to unrealistic fiber volume 
fractions. The effect of fibers on the macroscopic viscosity has to be factored 
into the macroscopic process simulation, as fibers are advected based on this 
pre-computed velocity field. Further, the regularization can lead to ’stuck’ 
fibers in regions with high shear rates [63]. 


2.2.6.2 Stokesian dynamics 


Most direct simulation approaches compute the motion of suspended parti- 
cles at low Reynolds numbers (Re « 1, Stokes flow) using an equation of the 


type 
My =f" +f" +f’ (2.65) 
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with particle mass M, particle velocity vp, hydrodynamic interaction force 
f, non-hydrodynamic interaction forces f! and stochastic forces causing 
Brownian motion f P [64]. The Brownian motion is neglected for fiber sus- 
pensions under investigation here and the hydrodynamic interaction is given 
as 

f? - - Ry (vp - v. + RD’ (2.66) 


with resistance tensors Ry, R; and the surrounding fluid velocity vo; [64]. 
Equivalent equations apply for the evolution of the rotational velocity of a sus- 
pended particle. In case of a suspended sphere with radius R, the resistance 
tensors take the form 

R,=6nnRI, R,=0. (2.67) 


Stokesian dynamics models do not account for two separate phases with an 
interface, but use the resistance tensors to capture all effects at a length scale 
smaller than the suspended particle. In most cases, the motion of fibers is 
driven by discrete evaluation of a given flow field væ, which is undisturbed 
by the presence of fibers (i.e. no long-range interactions). 


Yamamoto and Matsuoka [65] developed a beadchain model, in which a flexi- 
ble fiber is made up of several connected spheres (see Figure 2.7a). The fiber 
can bend, stretch and twist through linear elastic forces between spheres. 
Each sphere experiences a hydrodynamic drag according to Equations (2.66) 
and (2.67) as well as an equivalent torsion moment from a given fluid field [65]. 
They could demonstrate agreement of rigid fiber motion with Jeffery’s Equa- 
tion (2.21) and correct deformation of flexible fibers compared to Forgacs 
and Masons results [10]. 


The step from single fibers to fiber suspensions requires a model for short- 
range fiber-fiber interactions. A lubrication model for the normal direction 
between two rigid rods (see Figure 2.7b) was introduced by Yamane et al. [66] 
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in their simulation of a semi-dilute rigid rod suspension. They proposed the 
analytical solution based on a lubrication approximation 


2 |h 
pae e 


where | ha " describes the rate at which the cylinders approach each other 


at the contact area. They report a rather weak effect on fiber orientation 
(small Cj) and viscosity in the semi-dilute regime. Similarly, Yamamoto and 
Matsuoka extended their beadchain model to account for multiple fibers 
in a periodic cell with lubrication forces [67]. Lubrication forces only do 
not prevent fibers from penetrating each other at high aspect ratios. There- 
fore, mechanical contact forces have been added by Sundararajakumar and 
Koch [68]. The interaction of rigid rods and boundaries was investigated by 
Thomasset et al. [69]. 


For flexible fibers, the bead chain model was improved with the use of 
spheres [70,71] (see Figure 2.76), spheroids [72] (see Figure 2.7d), rods [73,74] 
2.7f) or combined beads [75, 76] (see Figure 2.7c) connected by sockets and 
joints to reduce the computational effort for fibers of long aspect ratios. Such 
models were applied to model flocculation [77], to obtain effective suspen- 
sion properties [78-80], analyze the effect of elastic fiber bending on the 
suspension elasticity [81], to model fiber jamming [82], to investigate fiber 
fracture [75] and fiber buckling [83]. 


Lindstróm and Uesaka [85] compute the velocity field on a grid instead of 
using a prescribed velocity field. The velocity field is solved including a body 
force field that opposes the hydrodynamic interaction forces (momentum 
conservation) and therefore includes long-range hydrodynamic interaction. 
They applied their two-way coupled approach to determine effective prop- 
erties of semi-dilute fiber suspensions [86, 87]. However, they utilize the 
drag of prolate spheroids, which leads to a total drag force on a fiber that 
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(a) Beadchain [65,67,84] (b) Rigid rods [66, 68,69] and (c) Combined beads [75,76] 
prolate spheroid [78,81] 


e nd 


(d) Socket and joint (e) Socket and joint (spheres) (£) Rod chain [73, 74] 
(spheroids) [72] [70, 77] 


Figure 2.7: Schematic overview on different approximations of fibers for Stokesian dynamics. 


depends on discretization [88]. These are the only models in this section, that 
compute a disturbance of the fluid field and therefore account for long-range 
interactions. 


More recent developments are targeting applications beyond representative 
volume elements further towards the component scale. For example, Kuhn et 
al. [89, 90] investigated rib filling during LFT compression molding. Hayashi 
et al. [91] suggested the use of constrained beams at this scale. Sasayama et 
al. [92] simplified the original beadchain model by dropping all rotational 
motion components and apply the model on the injection molding process 
of a center-gated disk with 20 wt% glass fibers [93]. 


2.2.6.3 Slender body theory 


Models using slender body theory are closely related to Stokesian dynamics. 
In this classification, such models solve an integral boundary equation for 
fibers without physical thickness and thus without a moment induced by 
the fluid flow. Further, slender body theory enables the computation of the 
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disturbance of the fluid flow by Green's functions of the Stokes equation 
called Stokeslets. 


For small Reynolds numbers Re « 1, an incompressible Newtonian fluid, 
vanishing velocity at infinite distance and a singular body force field, Equation 
(2.17) becomes 


0 = - grad(p) + ndiv(D) + ô (æ — xo) fo, (2.69) 


with a singular force fo in position x and the Dirac distribution 6(a — a9). 
The fundamental solution to this equation is given as 


1 I — £o) 8 (£ — 
T'(a- £0) = QUEE - a (2.70) 
8 | |æ- xo] |æ- 2o | 
and called the Oseen tensor. It is used to compute the Stokeslet 
v= fol (x — £o) (2.71) 


as a solution for velocity field [94,95]. An extension of this method for singular 
forces is obtained by distributing several Stokeslets along a line zo (s) with arc 
length coordinate s. The solution is then given by the integral 


J; 
v=f fo(s)F (a — zo(s))ds (2.72) 
0 


with a force per unit length fo(s) and is called the single-layer formulation for 
Stokes flow [40,95]. The linear superposition of this formulation is commonly 
used to represent slender, one-dimensional bodies suspended in a highly 
viscous flow. 


Hinch developed an early model for the motion of a single perfectly flexi- 
ble, but inextensible thread utilizing slender-body theory [96]. His model 
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provided early insight in fiber dynamics, but was not able to meet the ex- 
perimental results by Forgacs and Mason [10] due to neglected width of the 
thread by slender-body theory. 


The effect of hydrodynamic long-range interactions between rigid rods on the 
fiber orientation evolution was studied by Rahnama et al. [97] and showed 
good agreement with the results by Stover et al. [98]. Slender body theory 
faithfully computes the effective velocity increase up to a certain limit of 
fiber volume fraction, where it underestimates the effect [42] - likely due 
to neglected short-range lubrication interaction. Fan et al. [99] combined 
the slender body framework by Rahnama et al. for long-range interactions 
with the model for short-range lubrication by Yamane et al. They were able 
to determine effective rheological properties and Folgar-Tucker interaction 
parameters on periodic domains with rigid suspended fibers. The same 
method was applied to obtain an empirical fit for the Folgar-Tucker parameter 
depending on aspect ratio and fiber volume fraction [100]. A recent review 
on more advanced modeling of flexible fibers suspensions with a focus on 
slender-body theory is given by du Roure et al. [101]. 


2.2.6.4 Micro-macro approaches for highly concentrated 
suspensions 


This class of models is applied for the homogenization of highly concentrated 
fiber suspensions, in which long-rang interactions and viscous matrix stress 
may be negligible compared to short-range interactions [8, 102]. Le Corre 
et al. solve the momentum balance in a complex rigid fiber bundle network 
with viscous tangential lubrication forces and lubrication moments. The 
lubrication forces are formulated as 


1, g [levala 
7p. «il Gh e^ e 
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where the parameter G is interpreted as effective sheared gap and da is the 
major diameter of the elliptical bundle cross section [8]. The effective sheared 
gap G is different to the physical gap g in direction hag, since it describes a 
thickness that is equivalent to the complex shear zone between bundles and 
does not describe the physical distance between the bundle surfaces. The 
authors employ an up-scaling scheme to compute the macroscopic prop- 
erties of the suspension from the mesoscale model [8, 103]. Their analysis 
also showed that the lubrication moment is of minor importance for most 
cases and confirmed the observation by Servais et al. [57] that friction is of 
minor importance in bundle suspensions. Latter was also verified experimen- 
tally [104]. 


2.2.6.5 Resolved methods 


Resolved methods are direct numerical simulations that explicitly resolve 
the particle-fluid interface and solve the equations of motion for these two 
separate phases. These methods may be further subdivided in immersed 
boundary methods and particle methods. 


Immersed boundary methods use an Eulerian background mesh and a La- 
grangian body, enabling the method to properly capture the two-way interac- 
tion between those two phases. Such methods have been applied in the field 
of bio-mechanics to model the whirling motion of flangella to study bacterial 
locomotion [105], to investigate the flow around phytoplankton employing 
local mesh refinements [106] and multiple wood pulp fibers [107]. However, 
immersed boundary methods do not scale very well considering the number 
of fibers and the requirement of a volume discretization that must be smaller 
than the fiber diameter [101]. 


Particle methods are mesh-less Lagrangian methods to solve partial differen- 
tial equations and are interesting for direct simulations of fiber suspensions, 
because the phases can be naturally modeled by different types of particles. 
Kromkamp et al. [108] investigated the dissipation of suspended cylinders 
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with the Lattice Boltzmann (LB) method. The LB method was also applied to 
single flexible fibers and good agreement was achieved for bending modes 
and orientation orbits of stiff fibers [109]. Duong et al. applied dissipative 
particle dynamics (DPD) to Boger fluids, but several parameters have to be 
calibrated for correct hydrodynamic forces on fibers [110]. A moving particle 
semi-implicit method was applied by Yashiro et al. [111, 112] to simulate the 
injection molding process of dilute short-fiber-reinforced composites with 
connected particles. They can qualitatively reproduce molding effects, but do 
not verify the method and show only results with small fiber volume fractions. 
Smoothed particle hydrodynamics (SPH) was used by Yang et al. to model 
the motion of a single flapping fiber [113]. They represent the fiber by so- 
called element bending groups, which ensure an elastic connection between 
individual SPH particles. SPH was also used to model injection molding with 
multiple fibers [114, 115], but like Yashiro et al., the reported fiber number 
density of this work was unrealistically small. A more detailed analysis of SPH 
for fiber suspensions showed that it can be used to determine Folgar-Tucker 
constants on periodic domains, if a novel correction term for fiber thickness 
is added [116]. 


2.2.7 Concluding remarks on suspension models 


Macroscopic models utilizing fiber orientation tensors are numerically effi- 
cient models to describe dilute and semi-dilute suspensions with a fine distri- 
bution of short fibers. Such a suspension is macroscopically anisotropic due 
to long-range fiber-matrix interactions and short-range fiber-fiber interac- 
tions such as lubrication and contacts. There are different parameterizations 
for this anisotropic material response and several models are proposed to de- 
scribe anisotropy based on effective material properties as well as the second 
and fourth order fiber orientation tensor. However, the application of macro- 
scopic models is more suitable for injection molding process simulations, in 
which fibers are short in comparison to any geometrical feature. 
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Direct simulation models renounce the use of fiber orientation tensors and 
model individual fiber instances instead. The simplest approach within this 
class of models are kinematic models, where fibers are advected directly with 
a prescribed fluid velocity in combination with a regularization scheme to 
maintain a constant fiber length. These models are computationally fast, but 
do not affect the fluid flow and do not account for interactions between fibers. 
Resolved computational models with fully modeled individual matrix phase 
and fiber phase offer the potential to completely describe a suspension by 
immersed boundary methods or particle methods. However, these models 
are computationally expensive and limited to small representative volume 
elements. Slender body theory offers a computationally efficient approach 
to solve fluid flow and fiber motion simultaneously and is commonly used 
in the bio-mechanical community. While the method is very suitable for 
representative volumes with prescribed far-field boundary conditions, it can 
not be applied easily to component scale compression molding simulations 
with a deforming cavity and contacts at the mold surface. Similarly, the pre- 
sented micro-macro up-scaling approaches are restricted to representative 
volume elements. The most common approach for fiber suspension mod- 
eling is Stokesian dynamics which solves the force balance for interlinked 
suspended bodies due to hydrodynamic forces and non-hydrodynamic in- 
teraction forces. Within this class of models, the work by Lindstróm and 
Uesaka [85-87] is particularly noteworthy, as it affects the fluid flow to be 
anisotropic through two-way coupling. 


2.3 Sheet Molding Compounds 


Sheet Molding Compound (SMC) is a discontinuous fiber reinforced polymer 
made from uncured thermoset sheets and can be molded to parts under 
pressure and temperature by compression molding technology. Parts from 
this material offer a superior strength to density ratio (compared to polymers 
without reinforcement) at low cost. Furthermore, corrosion resistance and 
high achievable surface quality make it attractive for outer vehicle parts, such 
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as truck body panels, automotive hoods and trunks. While continuously 
reinforced polymers have limited forming capabilities due to the inextensi- 
bility of fibers in the fabric, SMCs can flow in complex geometries during 
compression molding. The incorporation of complex geometrical features, 
carbon fibers, metallic inserts and continuous fiber patches enables the appli- 
cation in structural applications such as subfloor structures in the automotive 
industry. 


2.3.4 Production 


The production of SMC parts is subdivided in two major steps: the production 
of SMC sheets from raw materials and the subsequent compression molding 
of sheets to a part. The sheet production (see Figure 2.8a )starts with mixing 
of all ingredients (thermoset resin, filler, thickener, additives, catalyst, mold 
release agent, inhibitor) except for fibers to prepare the paste that later forms 
the matrix material of the composite. The paste is then filled to doctor boxes 
of an SMC sheet production line, which apply the paste on two different 
plastic foils that run continuously through the machine. The boxes have 
adjustable plates to ensure a specific paste film thickness on the foil. 


At the same time, multi-end fiber rovings are pulled from several bobbins 
to a cutting unit that continuously cuts the strands in 25 mm or 50 mm long 
pieces. These bundles fall down on the paste layer of one of the plastic foils 
in a randomly in-plane orientated manner. The second resin-coated foil is 
placed on top of the first foil to enclose the resin and fiber bundles between 
both foils. The compound is then pulled through several (heated) calendering 
rolls to impregnate the fiber bundles with the thermoset paste and compact 
the compound. 


Finally, the sheet is rolled to a take-up coil and stored for several days and 
the resin starts a maturing process. During this maturing phase, the viscosity 
increases several orders of magnitude (compare Figure 2.8b). 


34 


2.3 Sheet Molding Compounds 


Roving Cutting 10° 
m 


Paste 
application 


Compression molding 


Shear viscosity in Pa s 
- 
e 
w 


É 
application Calendering leere 
0 1 2 3 4 5 
Time (days) 
(a) Schematic SMC sheet production line. (b) Exemplary shear viscosity [117]. 


Figure 2.8: SMC sheet production and maturing phase. The viscosity plot is recorded at ashear 
rate ý = 1.0s7! [117]. 


The sheets are cut, stacked and placed into a hot mold for compression 
molding. The stack typically covers a fraction of 2596 to 10096 of the mold and 
the remaining cavity is filled by flow of SMC, as the press closes the upper 
mold stamp in the corresponding lower mold half. The compression follows a 
defined closing velocity profile until the compression force hits a threshold at 
which the press controller switches to force control. This compression force 
level is held for a few minutes until the resin is sufficiently cured. Finally, the 
press opens the mold and the cured part is released. 


Possible defects according to Mallick [118] are summarized in Figure 2.9. If 
air is introduced in the paste during mixing or compounding and cannot 
be removed by compression between calendering rolls or if air is entrapped 
between sheets during stacking, pores and pinholes at the surface may form. 
Entrapped air or gaseous products of the curing reaction, which are closed 
during the compression may also form blisters that pop up as soon as the 
mold is opened. Sink marks may occur in resin rich areas due to shrinkage of 
the resin during curing. Typically, SMC is highly filled with chalk to prevent 
such surface defects. 


Fiber bundles reorient during the mold filling process and this may lead to 
unwanted fiber miss-orientation. For example, fibers close to mold edges 
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Figure 2.9: Possible defects in parts molded from SMC according to Mallick [118]. 


have a tendency to be oriented parallel to the edge. Additionally, knit lines 
may form where SMC has to travel around a feature. The line at which two 
flow fronts of SMC meet is typically weak, as only a limited amount of fibers 
crosses this knit line. If the paste travels faster than the fiber bundle network 
during molding, fiber-matrix separation (FMS) occurs. This phenomenon 
happens often in small ribs, at the flow front and other small features, where 
fiber bundles may get stuck. Fiber-matrix separation results in a certain area 
of a part consisting of matrix material only, which has significantly worse 
properties compared to the composite. 


2.3.2 Deformation mechanisms 


The micro-structure of SMC sheets implies certain restrictions on the de- 
formation kinematics, which are significantly different to the deformation 
of an isotropic fluid or solid. The initial stack is characterized by a random 
in-plane fiber orientation and the thickness is typically small compared to the 
lateral dimensions of the stack and the fiber length. Marker and Ford [119] 
present one of the earliest investigations on the deformation mechanisms 
and heating of SMC sheets. They investigate large stacks of ten sheets thick- 
ness in which the interlocking behavior of individual sheets resists transverse 
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shear deformation and outer sheets close to the mold wall experience larger 
deformation. Barone and Caulk [120] use this observation to simplify the 
kinematics of thin stacks (one or two sheets) by proposing that no shear 
velocity gradient exists along the thickness direction of the stack. This as- 
sumption was declared unrealistic by Tucker and Folgar [121], who suggested 
that a no-slip condition at the surface must apply. Subsequently, Barone 
and Caulk [122] found experimental evidence for their proposed stack de- 
formation mechanism using a disk shaped initial stack with black and white 
sheets of SMC. They performed compression experiments at two mold closing 
speeds (1.75 mms"! and 10.0 mm s^!) with multiple stack sizes. Exemplary 
results for a 5-layer stack are shown in Figure 2.10a and 2.10b. Barone and 
Caulk concluded that the deformation mechanism is primarily a uniform 
extension that can be accompanied by slip between layers and mold walls for 
large stacks and slow deformation speeds. Contrary to the fountain-flow ob- 
served in thermoplastic compression molding [123], this deformation mode 
is termed plug-flow. Later, Barone and Caulk modeled the slip at the mold 
surface with a constant friction, Coulomb friction and hydrodynamic friction 
and concluded that the hydrodynamic friction yields most realistic results for 
the SMC, they investigated [124]. 


Later investigations of the stack deformation differentiate a squish, flow and 
boiling phase during molding [125]. The squish phase refers to the initial 
complex flow front deformation due to thermal influences and release of 
air entrapped in the stack. Exemplary flow front shapes are illustrated in 
Figure 2.10c for 10 mms"! and in Figure 2.10d for 2mm s™} mold closing 
speed. Especially in the case of faster compression, it becomes apparent that 
the bottom layer exhibits larger deformation upon start of the compression, 
as it is hotter due to a longer contact with the mold surface. The subsequent 
flow phase is characterized by a stable plug-flow and the final boiling phase 
refers to gas release after complete closing of the mold. A minimization of the 
squish phase yields to a more homogeneous flow and reduced the volume 
fraction of pores in the SMC part [126]. Fiber bundles in the SMC core stay 


37 


2 Fundamentals and state of research 


(a) Fast deformation of 5-layer SMC [122] (b) Slow deformation of 5-layer SMC [122] 
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(c) Fast squish deformation [125] (d) Slow squish deformation [125] 


Figure 2.10: Deformation mechanisms of SMC reported in literature. 
typically intact during molding and only fiber bundles in the lubrication layer 


near the mold surface disentangle during the flow [117, 127-129]. 


2.3.3 One-phase modeling approaches 
2.3.3.1 Generalized Hele-Shaw models 


Hele-Shaw models simplify the incompressible Stokes flow between two 
plates to a two-dimensional problem, if the lateral dimensions are large 
compared to the gap between theses plates. They can be applied to model 
the compression molding process between two flat plates with gap h, that 
close with velocity h. The pressure field in such a 2D-domain is given by 


div(Sgrad(p)) =h (2.74) 
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where the divergence and gradient operators refer to the 2D-domain and 


h/2 2 
S= [ “= E dz (2.75) 
h/2 


is a measure of mobility [130]. Here, the position zo refers to the position of 
the flow, where dv,/dz = dvy/dz = 0. The corresponding 2D-velocity field is 
obtained as 

v=- > grad(p). (2.76) 


For the special case of isotropic Newtonian material behavior and no-slip 
conditions at the mold walls, S = h?/12n applies and zy = 0, as depicted in 
Figure 2.11a. These assumptions were used in the earliest models of SMC 
compression molding [131]. However, the thermally induced viscosity change 
may lead to failure of the parallel squeezing assumption and more complex 
velocity profiles develop, as the viscosity near the mold surfaces decreases 
with increasing temperature, as shown in Figure 2.11b [132]. A dimensionless 
criterion to determine whether a parallel squeezing assumption is valid was 
developed by Lee et al. [133]. 
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(a) Hele-Shaw flow (b) Thermal effect in SMC (c) Plug-flow assumption 


Figure 2.11: Velocity distributions through thickness for different mold boundary conditions 
and thermal conditions. 


The Hele-Shaw model for compression molding was solved for arbitrary 
domains beyond simple analytical cases using the finite element method 


39 


2 Fundamentals and state of research 


[121]. The authors indicated the incorporation of the energy balance and 
Non-Newtonian viscosity in their simulation framework. 


Generalized Hele-Shaw models utilize the simplification for the flow between 
arbitrary curvilinear surfaces with a thin gap between them. They are not 
restricted to flat plates and the application to compression molding was 
shown by Lee et al. [130]. The proposed model included heat transfer and 
isotropic Non-Newtonian viscosity and was solved using finite elements. 
The comparison to experimental results showed a good agreement for thin 
charges, but differences for thick charges due to a neglection of extensional 
stresses. The model was later combined with tensorial orientation models to 
predict the orientation after compression molding [134]. For a Newtonian 
material with constant part thickness, the Generalized Hele-Shaw flow may 
be solved using the boundary element method, which is very efficient and 
does not require a complicated meshing procedure [135]. 


2.3.3.2 Plug-flow models 


The Generalized Hele-Shaw models employ a no-slip boundary condition 
at mold walls, which does not comply with the observed SMC deformation 
behavior of Barone and Caulk [122]. Thus, the boundary element methods 
were extended to account for the modified flow behavior depicted in Figure 
2.11c with slip and hydrodynamic friction at the walls that maintains a veloc- 
ity profile without bulk shear deformation [136, 137]. Later, non-isothermal 
modeling (similar to the illustration in Figure 2.11b) confirmed the validity of 
a flat velocity profile with a thin hot lubrication layer close to the walls [138]. 
This settled the dispute between Lee and Tucker (advocates of the no-slip 
model) and Barone and Caulk (advocates of the plug-flow model) in favor of 
the plug-flow assumption. 


As the simplest approach, the hydrodynamic friction can be described by 
a hydrodynamic friction coefficient A, i.e. the proportionality constant be- 
tween relative velocity and friction force. For sufficiently thin parts, that is 
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hn! (B2A) « 1, friction dominates the momentum equation. In this case, the 
mobility term in Equation (2.74) can be replaced with 


h2 


goe 
24 


(2.77) 


to use an equivalent form of the generalized Hele-Shaw model with a more 
accurate description of the flow kinematics [139-141]. Such a model was 
successfully applied to the simulation of complex automotive structures as 
early as 1990 using finite elements [142]. The boundary element formula- 
tion requires too many simplifications for the application to such complex 
problems. This approach was extended for non-isothermal conditions later 
on [143]. Kotsikos et al. [144] employed a variational approach to allow shear 
and extensional deformation components and compared the results to exper- 
imental data from squeeze flow compression tests. They confirmed that SMC 
flow is predominantly of an extensional nature. 


After the establishment of basic models for the compression molding sim- 
ulation, the research focus switched to characterization of the material pa- 
rameters in order to predict molding forces correctly. Rheometry of SMC is 
difficult as the material has to be considered inhomogeneous at the scale of 
typical rheometers (such as plate-plate rheometers). Hence, custom press 
rheometers were developed to compress a sufficiently large sample and iden- 
tify rheological parameters with underlying modeling assumptions. Le Corre 
et al. performed lubricated (i.e. a lubrication film was applied to the mold 
surface to eliminate friction between SMC and the molds) disk compres- 
sion experiments and shear experiments with temperatures, deformation 
rates and fiber volume fractions similar to industrial applications [145]. They 
showed that the SMC behaves as power-law fluid, if the underlying paste is a 
power-law fluid and that the material is highly anisotropic. They identify the 
need to model SMC as a transversely isotropic material. These results were 
corroborated by plane strain tests to measure the transverse force and a set 
of constitutive equations was suggested to describe non-linear anisotropic 
SMC viscosity [146]. First parameters for non-lubricated flow, especially the 
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hydrodynamic friction parameter, were obtained by a comparison of the com- 
pression force between different initial positions of the stack in rectilinear 
compression [147]. The hydrodynamic friction can also be characterized 
more accurately by employing pressure sensors along the flow path and it 
was shown that both, friction and bulk rheology, are relevant to model correct 
compression forces [148]. Spiral rheometers with pressure sensors along the 
flow path were suggested as compact alternative to rectilinear press rheome- 
ters [149]. However, they may be considered as a rather qualitative method to 
compare SMC processing capabilities than as an exact tool for characterizing 
friction parameters. 


Several authors argue that the SMC mold filling process can be considered 
isothermal except for the heating of the small lubrication layer at mold sur- 
faces [148, 150]. Shokrieh et al. [151, 152] disagree with this assumption and 
suggest a time dependent thickness of the lubrication layer that is controlled 
by the development of the temperature in the stack and requires no additional 
experimental parameters. 


An additional characterization technique is tensile testing. It enables the 
characterization of damage and breakage in cases where sheets of SMC are 
drawn to deep cavities instead of being driven by a squeezing motion between 
mold surfaces [153]. 


While most previous models assume SMC to be incompressible, 3D in-situ 
computer tomography showed that the closing of pores that originate from 
sheet manufacturing cause a macroscopic compressible behavior [154]. This 
effect is especially true for structural SMCs and Dumont’s equations [146, 148] 
have been modified to account for this compressibility [150, 155]. 


2.3.3.3 3D Models 
Most early simulation approaches simplify the SMC mold filling as a 2D pro- 
cess. This is appropriate for large shell-like structures, but has limitations, 


when it comes to ribs or other complex 3D features in SMC parts, which are 
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increasingly important for structural SMC applications. An initial attempt to 
model SMC mold filling with a 3D CFD model was done by Vahlund [156], 
even though the method was applied to a shell-like part. Kluge et al. [157] uti- 
lized a 3D simulation to accompany their rheological experiments and gain 
deeper understanding of the material response for different viscosity models. 
Motaghi et al. [129] used the commercial process simulation software Au- 
todesk Moldflow, and employed an RSC model for fiber orientation. However, 
they used an isotropic viscosity model and did not compare corresponding 
compression forces. Hohberg et al. [158] pointed out that boundary condi- 
tions used in Autodesk Moldflow at that time do not respect the plug-flow 
assumption and suggested an elongational viscosity model with appropriate 
boundary conditions in Simulia Abaqus/ Explicit. The model showed poten- 
tial, but suffered from numerical instabilities. As of 2021, the latest version 
of Autodesk Moldflow allows the application of wall slip models at all mold 
surfaces!, which makes Autodesk Moldflow a viable software to predict SMC 
compression molding at macroscopic scale. The ability to account for wall 
slip is also available in CoreTech System Moldex3D and has been investigated 
for thermoset injection molding [159]. 


2.3.3.4 Thermoviscoplastic models 


Initial models for the SMC response during disk compression also suggested 
the use of a visco-elasto-plastic model [160] instead of the predominantly 
purely viscous models suggested above. However, this initial work did not 
account for the plug-flow kinematics described later by Barone and Caulk 
[122]. Kim et al. [161] suggested a modeling approach similar to that of 
sheet metal forming describing the SMC as solid anisotropic material with 
Hill’s yield criterion. They characterized parameters from lubricated disk 
compression tests and validated their model by comparing a finite element 


1 The required improvement was specified by Sven Revfi and Nils Meyer. Sven Revfi suggested 
the improvement to Autodesk and the issue was resolved in the new release. 
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simulation with experimental data of filling a T-shaped rib. In contrast to 2D 
shell simulations, this approach modeled a 2D cross section of the rib and 
was able to predict the filling behavior. Later, the authors enhanced their 
model by incorporating heat transfer, a curing model and a simple approach 
to predict the resulting fiber volume fraction [162] and extended the model 
to 3D finite elements [163]. The numerical results showed good agreement to 
molding trials of T-shaped parts, that were molded with multiple colored SMC 
layers. The characterization of the parameters for the anisotropic flow rules 
is described by Lin et al. [164, 165]. In most SMC applications however, the 
elastic component is rather small and the plastic deformation has little effect 
on the material response. Therefore, the majority of SMC models assumes 
fluid-like behavior, as described in the previous sections. 


2.3.4 Two-phase models 


Most SMC models consider the material as a one-phase continuum, as de- 
scribed in the previous Section 2.3.3. However, segregation of fibers and 
matrix is a relevant phenomenon in SMC compression molding, particularly 
in advanced applications with ribs, knit line formation and long flow paths. It 
motivates the development of models that can account for different velocities 
of fibers and matrix material. Besides the direct numerical models described 
in Section 2.2.6, this has been addressed by macroscopic SMC models which 
are shortly described in this section. 


Dumont et al. suggested a superposition of two immiscible phases (fiber and 
matrix) occupying the same volume at isothermal and incompressible con- 
ditions [166]. They obtain two momentum balance equations from mixture 
theory with a viscous momentum exchange term and solve the equations 
on a shell finite element mesh. The momentum exchange term is chosen 
such that it reduces to the equations of flow through porous media, if fibers 
are held in place rigidly. A similar approach was developed by Perez et al. 
resulting in a general 3D Brinkmann model [167]. However, the weighting 
factor between the two extrema (flow through porous media and suspension 
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flow with equal velocities between fiber and fluid) had to be prescribed and 
a lubrication approximation (compare Figure 2.11a) was used. Later, they 
related the weighting factor to the fiber-fiber contact density to compute the 
fiber volume fraction evolution [168]. 


2.4 Unidirectional patch reinforcements for 
co-molding 


Unidirectional fiber reinforced composites offer superior weight specific 
mechanical performance in fiber direction compared to discontinuous FRPs 
and most metals. They are either manufactured by thermoplastic tape laying 
combined with thermoforming, by consolidation of non crimp fabrics that 
are pre-impregnated with a thermoset matrix (prepregs) or by infiltration of 
dry fabrics (RTM, WCM, VARI). 


However, these materials are rather expensive and the ability to form complex 
shapes is limited by the inextensibility of fibers. The combination of discon- 
tinuous SMC with local unidirectional reinforced patches in a co-molding pro- 
cess offers potential to benefit from advantages of each constituent [3]. Thus, 
this section focuses on unidirectional pre-impregnated thermoset patches 
that can be co-molded with SMC, as developed by David Biicheler in the GRK 
2078 [169, 170]. 


2.4.4 Production 


The production of UD-patches for co-molding with SMC is similar to the 
production of SMC sheets and can be performed on the same production 
line (see Figure 2.12a). Instead of a cutting unit, the fabric is placed in front of 
the compounding zone and is impregnated with the matrix material between 
top and bottom foils. For the application in co-molding, the material is 
exposed to a special heat treatment after the impregnation to rapidly develop 
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a highly viscous resin B-stage state in the material. The B-staged material 
has experienced polymer chain growth, but no cross linking of the thermoset. 
This treatment allows direct cutting and forming after the impregnation 
process and ensures enough stability of patches during molding [170]. 


Paste SMC stack 
application Tempering SSS 

UD S ' ? 
Fabric 


=> 


UN m. > 
: u 
Preforming 
Co-molding 
nn 
application Calendering 
(a) Schematic patch production line (b) Schematic co-molding process 


Figure 2.12: Patch production and co-molding process. 


Patches are stacked while still being tacky directly after impregnation. Sub- 
sequently, they are cut, preformed and heated. After 60 min at 80°C, the 
patches are in a stable shape close to its final geometry in the co-molded part 
and are released from the preforming mold. The UD-patches are molded 
together with SMC made from the same matrix resin in a regular compression 
molding mold (see Figure 2.12b). During compression molding, patches may 
be subjected to the defects illustrated in Figure 2.13. Displacement occurs 
if the SMC excites a net force on the patch that is larger than the friction 
between patch and mold. The SMC acts upon the patch with normal forces 
and a sticky tangential friction. Local deformation is caused by a combination 
of insufficient friction to hold the patch in place and a distributed load by 
the SMC causing permanent deformation of the patch. If this deformation 
exceeds the strength of the patch in a particular direction, fracture may occur. 
This occurs typically perpendicular to the fiber direction by matrix tension or 
shear, as fibers have a high tensile strength. 
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Figure 2.13: Possible defects in unidirectional patches during co-molding (top view). 


2.4.2 Deformation mechanisms 


The quasi-inextensibility of fibers and comparably low shear stiffness implies 
kinematic constraints on the deformation mechanisms of continuous fiber 
reinforced fabrics. These deformation mechanisms are usually classified 
as intra-ply mechanisms happening in a single patch layer and inter-ply 
mechanisms occurring between individual layers of patches. The intra-ply 
mechanisms can be further categorized in bending behavior, membrane 
behavior and in-plane shear behavior [171]. Unlike homogeneous isotropic 
materials, bending stiffness and membrane stiffness of engineering textiles 
may require a decoupled modeling approach, if the polymer matrix stiffness 
is extremely small and allows slip between individual fibers during bending 
load [171]. Inter-ply mechanisms describe slip between plies, slip between 
mold and plies as well as delamination. The friction between mold and patch 
has been previously investigated by Bücheler [169] for the material modeled 
in this work. 


The predominant deformation mode of biaxial fabrics is called pure shear or 
Trellis shear that is motivated by inextensible fibers in two transverse direc- 
tions with contacts acting as pivot points. In contrast, UD fabrics naturally 
deform under so-called simple shear assuming inextensibility in fiber di- 
rection and incompressibility. The simple shear deformation is depicted in 
Figure 2.14. 


The deformation gradient F' maps the initial configuration to the current 
configuration, as indicated in Figure 2.14 by black coordinate systems with 
and without superscript ’0’, respectively. The coordinate system, which is 
transformed by F is called the covariant reference system. However, the 
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ex 
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Initial configuration Current configuration 


Figure 2.14: Simple shear in a unidirectional patch. The covariant coordinate system is depicted 
in black (ex, ey) and the co-rotational system in green (éx, éy). 


deformation gradient includes a rigid body rotation and is therefore not 
objective. A polar decomposition 


F=RU=VR (2.78) 


in a pure rotation R and a right stretch tensor U or left stretch tensor V 
provides a separation in rigid body rotation and an objective strain mea- 
sure. However, the by R rotated co-rotational frame is typically not suited to 
describe the deformation of anisotropic composites, because the principal 
material orientations become non-orthogonal during shearing, as shown 
with the black coordinate system in Figure 2.14 in comparison to the orthog- 
onal frame depicted as green coordinate system. A possible solution to cope 
with this is the introduction of additional transformations in a hypoelastic 
constitutive model (see Section 2.4.3.1) or the formulation of a hyperelastic 
constitutive model (see Section 2.4.3.2) in the initial configuration. In any 
way, the constitutive model needs a deformation measure that does not de- 
pend on rigid body motion and rigid body rotation (objective strain measure). 
The definition of the right Cauchy-Green tensor 


C=F'.F=U? (2.79) 
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allows the definition of the Green-Lagrange strain 


E= C-D, (2.80) 


which is a suitable objective deformation measure for modeling the patches. 


2.4.3 Modeling approaches 


The simplest approach to model patch deformation are kinematic models 
that assume incompressibility and inextensibility of fibers [172]. In biaxial 
fabrics the contacts between weft and warp fibers act as pivot points with- 
out slippage and in unidirectional reinforcements ideal simple shear with 
constant fiber spacing is assumed. The forming result is purely controlled 
by geometrical considerations and without solving a momentum-balance 
(similar to the kinematic fiber models presented in Section 2.2.6.1). The 
benefit of such models is the fast computation time, which can be lever- 
aged in optimizations [173], but the prediction capability of defects is very 
limited [174]. 


A more rigorous approach is the constitutive modeling of fabric deforma- 
tion that accounts for material behavior. Such models are developed for the 
microscale to analyze slippage of fibers and resulting yarn cross sections 
during weaving [175]. Detailed mesoscale models [176, 177] and discrete 
models build from 1D or 2D elements representing constituents of the com- 
posite [178-180] are used to determine effective properties based on the 
underlying yarn structure or directly model deformation at component scale. 
Most commonly, forming behavior of engineering textiles is described by 
homogeneous macroscale models. This is often done either by formulations 
that relate the rate of the Cauchy stress to the strain rate (hypoelastic) or by 
formulating the material behavior with respect to the initial configuration 
(hyperelastic) [181, 182]. 
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2.4.3.1 Hypoelastic models in co-rotational frame 


Hypoelastic models follow the general form 


a’ =C:D (2.81) 


V such 


with a tangent stiffness tensor C and an objective stress derivative o 
as the Jaumann rate or Green-Naghdi rate [183]. The constitutive equations 
are formulated in terms of the current Cauchy stress, which is integrated in 
time by an appropriate numerical scheme [184]. In general, fiber reinforce- 
ment fabrics are described in a non-orthogonal frame and a transformation 
between the orthogonal objective Cauchy rate (see 65, êy in Figure 2.14 for 
Green-Naghdi frame) and the non-orthogonal fiber frame (see ex, ey in Figure 


2.14) has to be applied [185-187]. 


2.4.3.2 Hyperelastic models in initial frame 


Hyperelastic models allow for a formulation of constitutive material models 
with respect to the initial configuration and thus omit the need for difficult 
coordinate transformations. Such models propose the existence of a strain 
energy potential W and compute the second Piola-Kirchhoff stress as 


ðW (C) 
9C ' 


S=2 (2.82) 


where C is the right Cauchy-Green tensor [183]. The strain energy potential is 
typically formulated in terms of invariants of the right Cauchy-Green tensor 
and includes energy stored by tension and shear [188, 189] as well as coupling 
terms [190,191]. 


The Cauchy stress can be computed by a push forward transformation with 
the deformation gradient as 


o= SF-S-F", (2.83) 
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where the Jacobian determinant J = det(F') denotes the volume change [183]. 


Given the non-orthogonal deformation behavior illustrated in Figure 2.14, a 
hyperelastic formulation w.r.t the initial configuration seems most appropri- 
ate to model the deformation of unidirectional patches during co-molding. 
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3.1 Research hypothesis 


Fibers in SMC are typically 25 mm to 50mm long and are required to fill 
geometric features of similar dimensions in advanced (semi-)structural appli- 
cations. This violates several conditions for the application of tensor based 
fiber orientation descriptors, e.g. fibers are not straight, fibers are not much 
shorter than geometric features and the effects are not local, because the 
motion of one fiber end affects the other end up to 50 mm away. The question 
is in which scenarios these violations may still lead to reasonable results and 
which alternative model could be used to alleviate these issues while still 
being applicable at component scale with reasonable computational times. 


The main hypothesis is that a mesoscale process model that accounts for 
the motion of individual fiber bundles improves the prediction quality of 
fiber architecture in SMC compression molding simulations at component 
scale. The improved prediction of the architecture can be qualitatively and 
quantitatively validated by comparison to uCT scans. 


3.2 Objective 


The objective of this work is the development of a 3D numerical simulation 
method for the compression molding process of SMC in regions, where fiber 
orientation tensors may not be suitable descriptors of the fiber orientation 
state. These regions are geometrical features such as narrow ribs and beads 
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which are common in advanced structural SMC applications as well as the 


area near mold walls. The goal is the prediction of structural defects during 


SMC compression molding, such as knit line formation, fiber-matrix separa- 


tion and fiber misalignment. Additionally, accurate predictions of necessary 


compression forces and fiber orientation states are desired. Further, the sim- 


ulation method should be able to incorporate a model for continuous fiber 


patches to predict processing effects in parts made from CoDiCoFRP. 


The process step under consideration is the mold filling process starting from 


positioning of the initial SMC stack and patches until complete fill of the 


mold. Preforming of the material placed in the mold is not considered, nor is 


subsequent curing, cooling and warpage of the part. 


The detailed sub-objectives are summarized as follows: 
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Development of a macroscopic SMC reference model utilizing a tensor 
based description of fiber orientation states to be able to compare 
mesoscale simulations to corresponding homogeneous macro models. 


Development of a mesoscale SMC model based on direct numerical 
simulation of the bundle suspension. The model should account for 
two-way bundle-matrix interaction, bundle-bundle interaction and 
confinement by molds. 


Development of a suitable mechanical model for unidirectional rein- 
forcement patches. 


Experimental characterization of required parameters for the devel- 
oped models. 


Validation of models by comparison to experimental process data and 
CT analysis of the resulting parts. 


Integration of the models in a CAE chain workflow to be able to transfer 
results from the process simulations to subsequent warpage or struc- 
tural simulations. 


3.3 Outline 


3.3 Outline 


The work is laid out as follows: The next Chapter 4 formulates a 3D macro- 
scopic reference model for SMC employing fiber orientation tensors and 
anisotropic viscosity. It is solved using the Coupled Eulerian Lagrangian 
framework in Simulia Abaqus/Explicit and a custom VUMAT subroutine. The 
equations are reduced to the 1D case in a press rheometer to investigate 
principal effects of individual parameters and provide a tool for rapid compu- 
tation of reference results. 


Chapter 5 includes the main contribution of this work, which is the develop- 
ment of a mesoscale model that is capable of running at component scale to 
provide more detailed prediction of defects and features during compression 
molding. The chapter is subdivided in the derivation of long-range bundle- 
matrix interactions and short-range bundle-bundle interactions. Sub-models 
are used in each case to derive phenomenological expressions for relevant 
terms at mesoscale, such as the drag coefficients of bundles and the effective 
sheared gap between bundles. Finally, a verification ensures that the model 
is able to reproduce analytical reference solutions. 


A model for unidirectional reinforcement patches is described in Chapter 6. 
The model describes each patch ply with two shell layers with shared nodes. 
One layer represents the unidirectional bundles and the matrix, the other 
layer models the unidirectional stitching yarn. Both layers are hyperviscoelas- 
tic and a Hashin damage criterion is employed to describe failure of patches 
during molding. 


Chapter 7 describes the characterization of necessary parameters for the 
SMC model and the patch model in the uncured state. This includes the 
transverse heat conductivity of SMC sheets and patch plies as well as the heat 
conductance from a steel mold to the SMC stack and patch. Further, the rate 
dependent and temperature dependent viscosity of the SMC paste and the 
compressibility of a uncured stack of SMC are determined. For patches, the 
relevant anisotropic visco-elastic properties are determined. 
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In Chapter 8, several experimental results and corresponding simulation 
results are compared to analyze the merits of direct mesoscale simulations 
and macroscale models. The first case investigates the general ability of the 
DBS to predict fiber architectures under quasi-incompressible isothermal 
Newtonian conditions in a double curved dome. The resulting fiber orienta- 
tion, curvature and knit line formation agree well with computer tomography 
scans of experimental specimens. The second case investigates fiber-matrix 
separation at ribs with quasi-incompressible isothermal Newtonian model- 
ing conditions and compares the results to pyrolysis results of specimens. 
The full non-isothermal non-Newtonian compressible model is applied to a 
press rheometer simulation in the next example. The resulting compression 
force and pressure at sensor positions is compared to experimental data to 
evaluate the model’s capability to not only predict fiber architectures, but 
also to predict realistic compression forces. The full model is also applied to 
a more complex small demonstrator part featuring ribs, beads and unidirec- 
tional reinforcement patches to show the ability to integrate the patch model 
with DBS. The final validation example is a large part to demonstrate that the 
suggested approach can be applied at the scale of larger SMC components. 


The modeling domains and investigated features are visually summarized in 
Figure 3.1. 
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Patches: S3RT/S4RT Mold: R3D3/R3D4 
- Hyperviscoelastic model - Hydrodynamic friction to SMC paste 


- Simple damage model - Virtual press controller 


- Non-isothermal - Rigid and isothermal 


Chapter 4 
Paste: EC3D8RT Bundles: T3D2 
- Non-Newtonian - Flexible/thread-like 
- Compressible - Realistic volume fraction 


- Non-isothermal - Drag from paste motion 


- Body force field from bundles - Interaction between bundles 


Chapter 4/5 


Figure 3.1: Graphical overview of the modeling domains and corresponding simulation fea- 
tures. The molds are represented by rigid shells (R3D3/R3D4) at a homogeneous isothermal 
temperature. The motion is controlled by a virtual press controller that mimics the controller 
of the physical press with a switch from displacement control to force control. A key feature is 
the explicit model of fiber bundles by flexible thread-like truss elements (T3D2) at component 
scale and with realistic fiber volume fractions. The fiber motion is governed by drag forces from 
the fluid phase and fiber bundles interact with each other as well as with the mold surfaces. 
The fluid phase is modeled in an Eulerian frame (EC3D8RT) and is generally non-Newtonian, 
non-isothermal and compressible. The drag of fiber bundles is subjecting the fluid phase to a 
body force field that causes anisotropic flow behavior. The simulation framework may incor- 
porate continuously reinforced patches modeled as shells (S3RT/S4RT), which are modeled as 
hyperviscoelastic material with a Hashin criterion for damage initialization. 
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4.1 General three-dimensional model 


Typically, SMC process simulation is conducted at a macroscopic continuum 
scale, which efficiently describes the flow of the SMC fiber suspension as 
a homogeneous material with effective material properties. To verify the 
mesoscale model developed subsequently in Chapter 5, such a macroscopic 
model is formulated here. It is implemented in the commercial finite element 
software Simulia Abaqus/Explicit and provides a reference solution to discuss 
differences in the modeling scales. 


4.1.1 Assumptions and scope 


The homogeneous description of SMC implies some assumptions, which are 
introduced to simplify the chemo-thermo-elasto-visco-plastic behavior of 
SMC: 


* Scale separation applies and the fiber architecture is sufficiently de- 
scribed by the second order fiber orientation tensor 


e SMC is a homogeneous one-phase material with effective properties 
* Curing does not affect the flow process 


* The deviatoric mechanical response of the effective material is purely 
viscous 
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The first two assumptions are reasonable for simple plate-like structures 
with in-plane dimensions that are much larger than fiber length. However, 
these assumptions are highly questionable in confined regions, knit lines 
and any feature at fiber length scale. The third assumption is founded in the 
hypothesis that flow occurs on a time scale much faster than curing (10° s vs. 
10? s). Thus, curing is neglected to reduce the number of parameters in the 
model and characterization efforts. However, even a small onset of curing 
might affect the viscosity significantly and thus introduce errors in the simula- 
tion results. The fourth assumption reduces the thermo-elasto-visco-plastic 
material behavior to a thermo-viscous behavior. There is certainly some 
elasticity in the composite (fibers are elastic and the polymer matrix induces 
entropy-elasticity as well as some early degree of cross-linking). Even so, the 
elastically stored energy in the SMC is assumed to be small compared to the 
dissipated viscous energy during the molding process. Not all dissipated 
energy is necessarily rate dependent. Internal friction, such as between fibers, 
causes some degree of plastic behavior. Contrary to fiber suspensions with 
dispersed bundles, which show a Herschel-Bulkley behavior with yield stress, 
fiber bundle suspensions are dominated by viscous fiber interactions [56, 57]. 
Hence, plasticity is neglected here for simplicity. 


4.1.2 Governing equations 


The initial boundary value problem is formulated on a domain x € D c R?. 
The solution fields are the velocity v : D — IR?, the mass density p : D — 
R, the temperature T : D — R and the inner variables describing the fiber 
orientation state A : D — A. Considering symmetries, this leads to a total of 
eleven independent variables (vi, v2, v3, p, T, A11, A22, A33, A12, A23; A13}. 


The solution is obtained from the macroscopic mass balance 


dp 


— - di - 0, 4.1 
ar iv(pv) (4.1) 


60 


4.1 General three-dimensional model 


the macroscopic momentum balance 


g 2 + div ((pv) e v) = div(o), (4.2) 


the macroscopic balance of inner energy 


O(pcp T) 


ao div((poy T) -v) = -div(d) +0 : D, (4.3) 


and an evolution equation for the second order fiber orientation tensor 
ôA 
5, divAev)- W-A-A-W+A(D-A+A-D-2A:D). (44 


The employed orientation model is the simplest possible version (see Section 
2.2.4) without diffusive interaction parameters. This could be easily replaced 
by a more complex orientation model, if the required additional parameters 
are available. However, it is shown later in this work that the simple base 
model shows remarkable agreement with the detailed mesoscale model for 
planar SMC compression molding. 


The governing equations result in a total number of eleven equations for the 
unknowns. The heat flux d and momentum flux ø are related to the solution 
variables via the constitutive equations outlined in the next section to close 
the problem. 


4.1.3 Constitutive equations 


The Cauchy stress is computed according to Equation (2.33) from several 
individual contributions as 


o —- -p(^e) +04 + ony t+ Opp (4.5) 
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where the pressure p is no individual solution variable, but modeled with 
an equation of state that depends on the current mass density through the 
Hencky strain 


se=tog( 7} = 1022] (4.6) 
Vo p 


A relation between mass density and pressure is determined in Chapter 7. 
Alternatively, a sufficiently large bulk modulus may be employed to model 
quasi-incompressible SMC. The isotropic deviatoric matrix stress 


cw -729(05 y): D' (4.7) 


uses a matrix viscosity that generally depends on shear rate y and temperature 
T. The extra stress from long-range fiber-matrix interactions is 


Afrin(T y) 
Oo = 
FM 3[In@/f)+InInd/f)+ C] 


; Ie A :D' 
(4.8) 


2 
+2nsn(T, y) | AUJ «T A-SIeAJ:D' 


as given in Equation (2.39) in combination with the parameter 7», sr from 
Equation (2.42). Strictly speaking, the resulting Equation (4.8) is only valid 
for Newtonian fluids with constant viscosity due to the assumptions made 
in the underlying derivations. It may be applicable, if the non-linearity in- 
troduced by (7, y) is not too severe (approximately constant local viscosity 
at the length scale of fibers), but results have to be reviewed critically when 
using Non-Newtonian viscosity. The shear factor ns is included to improve 
numerical stability in 3D models and might be considered a numerical pa- 
rameter that is sufficiently small to not affect the result but large enough to 
suppress spurious shear. The last term in Equation (4.5) was obtained by pro- 
jecting Equation (2.59) with the deviatoric projector P, to ensure that it does 
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not generate volumetric stresses. Hence, this extra stress from short-range 
fiber-fiber interactions is formulated as 


Ak g 1 
Bap = Pere) | — „IeB] :D'. (4.9) 


The thermal flux for the energy balance in Equation (4.3) is modeled accord- 
ing to Fourier’s law 
d=-xgrad(T) (4.10) 


where & represents the heat conductivity tensor. The conductivity is reduced 
to a scalar value « = «I here, to reduce the number of parameters to be char- 
acterized. This scalar model is a simplification, as in-plane heat conductivity 
is different to the transversal conductivity. However, the conduction process 
in SMC is controlled by the heating between hot molds transverse to the 
sheets, while the majority of in-plane heat transport is of convective nature 
due to the material flow. 


Finally, the fourth order orientation tensor À is computed from the IBOF 


closure [21] and the fourth order interaction tensor B is computed from the 


quadratic closure provided in Equation (2.63). A summary of the parameters 
for the macroscopic model is given in Table 4.1. 


Table 4.1: Parameters for the macroscopic SMC model. Additional parameters are needed to 
describe the relations for n(T,y) and p(p). 


Description Symbol 
Fiber volume fraction f 

Fiber aspect ratio Tp 
Shear viscosity of the matrix n(T,y) 
SMC compressibility relation p(Ae) 
Short-range interaction parameter kp 
Shear factor Ns 
Transverse heat conductivity K 
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4.1.4 Friction model at mold surface 


As mentioned in Section 2.3.3, Barone and Caulk investigated Coulomb fric- 
tion, constant friction and hydrodynamic friction as candidates for the thin 
lubrication layer between SMC core and the mold surface [124]. Hydrody- 
namic friction became the most common model and is generally described 


by 
m-1 
m=-a| l] Us, (4.11) 
Uo 


where 7y is the stress at the mold surface, A is a hydrodynamic friction 
coefficient, m is a power-law coefficient, vo is an arbitrary reference velocity 
for non-dimensionalization of the power-law term and vs is the slip velocity 
at the mold surface [148, 150]. The determination of friction parameters in 
literature focused mostly on a flow phase, where a stable plug-flow has formed 
and the initial squish flow has not be investigated in the same depth. During 
the squish phase, hydrodynamic friction models have the disadvantage to 
prescribe no friction stress at zero velocity. However, the initial contact may 
berather described by a sticking behavior in some SMCs, which would favor a 
model with constant friction stress To. Such a model restricts the slip velocity 
to v, = 0 for wall shear stresses smaller than ro. 


In this work, a transition function T is suggested to smoothly transition be- 
tween a constant value ro, which is applicable in absence of relative motion 
and a hydrodynamic friction model for relative motion. The transition func- 
tion is defined as 


r- ol x (4.12) 


with a parameter v; that describes the transition velocity between sticking 
and slipping states. The maximum shear stress at the mold surface is then 


m-1 
m=-rrfo] -a-mat Us. (4.13) 
0 


formulated as 
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Constant friction 


— Linear hydrodynamic friction 
Power-Law hydrodynamic friction 
Transition 


Mold wall stress || 74 ll 


Relative velocity at mold surface || vs || 


Figure 4.1: Comparison of friction models depending on relative velocity. While hydrodynamic 
models are widely used in literature [124, 148, 150], the initial squish phase may require a sticking 
constant friction at low slip velocities. Therefore, a transition model for sticking SMC is proposed 
in this work. 


A graphic summary of constant friction (v; — oo), linear hydrodynamic fric- 
tion (m = 1, v, — 0), power-law hydrodynamic friction (m # 1, v; —^ 0) anda 
transition model is presented in Figure 4.1. 


4.1.5 Virtual press controller 


The physical press follows a press profile given as M pairs of mold gap and 
corresponding closing velocity, as shown in Figure 4.2. Eventually the press 
controller switches to force-control in order to limit stresses on mold and 
press. A virtual press-controller is used to mimic this behavior as boundary 
condition during the compression molding simulation. 


As long as the compression force is below the force at switch-over Fmax, the 
profile is linearly interpolated to obtain the current press velocity. Regions 
outside the prescribed profile are linearly extrapolated. The linear interpo- 
lation/extrapolation w.r.t the gap leads to a non-linear velocity profile w.r.t 
compression time. 
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Figure 4.2: The press profile is given as M pairs of mold gap and press speed and is linearly 
interpolated w.r.t the gap. 


After the switch-over at time increment n,, a simple discrete PI-controller is 
employed to determine the current press speed 
n 


AE 
fins = hn + Ppen Pi >, —At (4.14) 


q-ns 
from the normalized errors 


Fmax- Fa; 
7 HM» 


AEg= (4.15) 


Fmax 


where hy describes the final compression speed of the press. The normal- 
ization ensures reliable force-control through a wide range of simulation 
parameters with constant control parameters Pp = P; = 0.5. 


4.1.6 Solution procedure 
The Coupled-Eulerian-Lagrangian framework of the commercial Finite Ele- 


ment solver Simulia Abaqus/Explicit is employed to solve the equations, as 
the deformations during compression molding are too large to be described 
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in a Lagrangian frame. This operator splitting method divides the equations 
in a Lagrangian part and an Eulerian part [192, 193]. The Lagrangian part 


SE sii 4.16 
xL (4.16) 
Ope). dite) (4.17) 

at |; 

ð T 

Wael) 2er (4.18) 

0t | 
u =W.A-A-W+A(D-A+A-D-2A:D) (4.19) 

1 


contains only the source terms and neglects the convective terms. These 
equations are solved in a traditional Lagrangian frame by explicit time inte- 
gration and enable to account for contact with other Lagrangian bodies. A 
subsequent step is used to solve the convection problem 


OP) a divans (4.20) 

t |E 
g div((ev) ® v) =0 (4.21) 

at |p 

O(pcp T) u 

3 ,rdiv(e mv) -0 (4.22) 
—| +div(Aev) - 0. (4.23) 

E 


This step effectively moves back nodes to remain an Eulerian frame and 
computes the flux through element boundaries by a second order advection 
scheme [192, 193]. An illustration of the process is given in Figure 4.3, where 
a circular shape is stretched horizontally. The initial shape is represented by 
a volume fraction of SMC e € {R |0 < e < 1} in each element and properties 
of partially filled elements are volume averaged. 


The contact between Lagrangian bodies and the Eulerian material surface uti- 
lizes a material reconstruction algorithm [194]. The surface is approximated 
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Langragian step wA 


Figure 4.3: Illustration of the CEL steps. The first step is equivalent to a conventional Lagrangian 
step, if the spatial time derivative would be replaced by a material time derivative. The second 
step moves the mesh back to its original configuration, computes the flux of material through 
element faces and adjusts the Lagrangian variables [193]. 


by a set of planes for each element that does not necessarily form a continu- 
ous surface. The reconstruction allows a definition of contact overclosures 
for a set of seed points on the contacting surfaces and thus allows to compute 
contact forces from a penalty approach. The resulting local contact force 
is then distributed on the elements nodes with their corresponding shape 
functions [193]. 


The constitutive model is implemented in a VUMAT subroutine, the virtual 
press controller is implemented as VUAMP subroutine and the friction is 
modeled with a pre-computed look-up-table depending on the absolute 
value of slip |v. |. 


4.2 Reduction to one-dimensional plug-flow 
in a press rheometer 


In some cases, such as the elongational plug-flow in a one-dimensional press 
rheometer, the previously described equations may be further simplified to a 
one-dimensional set of equations. It is desirable to be able to compute fast 
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solutions to such a problem to verify more complicated three-dimensional 
implementations and to identify parameters with a press rheometer. There- 
fore, this section states a set of one-dimensional PDEs to approximate the 
compressible, non-Newtonian plug-flow in a press rheometer. Subsequently, 
an analytical solution for the further simplified incompressible Newtonian 
case is derived. 


Figure 4.4: Dimensions of a press rheometer for one-dimensional elongational flow in the initial 
configuration (t = 0). 


The dimensions of a generic press-rheometer are sketched in Figure 4.4. The 
solution domain is then described by a non-dimensional scalar position x* — 
x/ X, where the position x is normalized with the flow front position X and 
x*eE {R |0<x*< 1}. The fiber orientation is assumed to be planar, such that 
only the Axx, Axy, Ayx and Ayy components of A remain. Further, all entries 
of A except for Ayxxx, Ayyxx, Axyxx and corresponding symmetries vanish. The 
vertical strain rate is prescribed by the deformation speed D,, = h/h in terms 
of the mold gap h. 


Isothermal bulk flow is a common assumption in literature (see Section 2.3). 
However, even small changes in temperature can have significant effect on 
the matrix viscosity. Hence, an analytical solution is employed to compute an 
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approximate average mold temperature during compression. The tempera- 
ture between two closing plates with constant closing speeds and ideal heat 
transfer at the mold surfaces (Dirichlet boundaries) can be expressed as 


T(t, Z) = Ty + (Tr — To) 


" — (4242 
2 cos(qz — 1) pega). | M (4.24) 


T 2 q h hohpcp 


with a constant mold surface temperature Tr, initial homogeneous stack 
temperature Tọ and q € N [132, 140]. Averaging the temperature distribution 
over the thickness 


h 
T(t)- Ji T(t,z)dz (4.25) 
h Jo 


yields 


T(t) = Tr + (Tr — To) 


oo _ coda 
1 > cos(qz — 1) : | qu <) (4.26) 
T 

q=1 


e 
q? hohpoy 


as an approximation of the average temperature in the SMC stack. 


4.2.1 Scalar set of PDEs 


For a model without short-range interactions, the constitutive relation may 
be formulated as 
cg =-p(p)I+V:D' (4.27) 


with an anisotropic viscosity tensor 
=. =. 1 
V=2n(T,y)P2+ n», se CT, y) [A-5184] (4.28) 


from equations (2.35), (2.39) and (2.42). With the assumptions stated above, 
this reduces to 


4 - = 1 
Va = 30053) + Ne, se 5) [aa - An (4.29) 
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2r : 1 

Vazz = E Y)+n2,sr(T,y) {Anca sa) (4.30) 
2 = - la 

Vozxx = -an bye ga, se (T, y) {Ara 3 Au (4.31) 
Ay oe = 1 

Vzzzz = ich y)+n2,sr(T,y) [Aree 73 (4.32) 


for the relevant terms. The crossed out components vanish due to the as- 
sumption that all fibers are in the xy-plane. 


Mass balance, momentum balance and the orientation equations are then 
formulated on the normalized domain as 


6 = -pDx- pDzz (4.33) 
ð T 

pv= Ax* (—p(p) + Vox Dyx + Vsxzz Dzz) E 27 (4.34) 

Axx = 2(Axx — Axxx) Dxx (4.35) 

Ayy = —2AyyxxDax (4.36) 

Axy = (Axy — 2 Axyxx) Dax» (4.37) 


where domain remains undeformed are progression of the flow front is 
OC) — 1 de) 

CX Ox’ 
the momentum equat describes a friction stress t that acts on both mold 


account for by 5 


e.g. Dy = fe x. The additional source term in 
surfaces (hence the factor 2) against the flow direction. The unknown fields 


are the mass density p, velocity v and the orientation components Axx, Axy 
and Ayy. 


4.2.2 Initial conditions 
The stack is considered to be initially at rest and at rest density po 


Y\,-9=% Plı-o= Po (4.38) 
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because it lays unconstrained on the lower mold half. The orientation is 
considered planar isotropic 


Aylı 9 =05, Ajyl 9 0-5; [4 = 0.0, (4.39) 


if the SMC sheet production process (see Section 2.3) does not induce any 
preferential fiber orientation. 


4.2.3 Boundary conditions 


The Dirichlet boundary condition at x* = 0 is 


Ul ag =O (4.40) 


x*= 


at all times. The boundary condition at the flow front x* = 1 is given by the 


flux boundary 
—P(0) + Vx Dxx + VxxzzDzz eg = 0 (4.41) 
during flow and by the Dirichlet boundary 
vl, 70 (4.42) 


as soon as the mold is completely filled. All other boundaries are zero-flux 
boundaries. 


4.2.4 Solution procedure 
The set of coupled PDEs, initial conditions and boundary conditions is com- 
pleted with algebraic expressions for utility variables. The first expression is 


the evolution of the flow front position 


X=vle-ı (4.43) 
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for X < Xmax. The second expression is the compression force 
1 
F- Bf — P(e) + VizxxDxx + Vzzzz Dzz dx" , (4.44) 


which is integrated from the normal stress at the mold surface. 


The set of coupled PDEs is solved using two pdepe solvers in MathWorks 
Matlab. The first solver integrates the equations until the mold is completely 
filled and a second solver restarts at that point with modified boundary 
conditions at x* = 1. The algebraic expression for the flow front is solved 
simultaneously with a forward Euler time discretization and the compression 
force is evaluated using trapezoidal integration. The computation of the 
compression force is needed at run time, as it is used by a virtual press 
controller to switch from displacement control to force control, as soon as 
the maximum compression force is reached. 


4.2.5 Analytical solution 


Further simplification of the one-dimensional flow allows to derive a simple 
analytical solution. Assuming incompressibility (Dz; = - Dy, = D) leads to a 
direct solution for the flow front position 


ho 
Xe (4.45) 
h 
and the velocity field 
v=-Dx. (4.46) 


If the fiber orientation state is approximately constant, i.e. Aj = 0.375 and 
Axx = Ayy = 0.5, the set of equations simplifies to just the momentum balance 


ap D?x+2— (4.47) 
ax? h i 
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Integrating this equation once for Newtonian viscosity and a constant friction 
To yields the pressure field 


T 1 
p--pD*o?- X^) +2 = (x- X)-2mD- N2, sP (Au zAw)D, (4.48) 
N Nm nn nn 
from boundary condition 
where the last term originates from the boundary condition (4.41). Finally, 
the compression force can be computed from Equation (4.44) as 


2 2,3_ 10 y2 

Fe = B| Z pD’ X°— X? —AnDX ~ 12, s DX Axox (4.49) 
Nee es es S ——M 
inertia friction matrix fiber-matrix 


with a contribution from inertia, friction, matrix viscosity and extra viscos- 
ity due to fiber-matrix long-range interactions. The deformation rate D is 
typically negative (compression) and thus each component adds a positive 
contribution to the overall compression force. 


In case of a linear hydrodynamic friction r = Av, the solution is 


2 A 
F,=B 5PD’ X? DX 4nDX — 9o, sr DX Ax |, (4.50) 


which reduces to the result obtained by Castro [139, 140], if the inertial term 
and the fiber-matrix term are dropped. 


4.3 Verification of one-dimensional plug-flow 


An exemplary solution of the 1D PDE model is given in Figure 4.5. In this 
generic example, the mold is closed at a rate of 1 mm s^! from an initial gap of 
9 mm to a final gap of 2 mm. Due to compressibility of the material, the flow 
front advances in a non-linear fashion until the mold is completely filled after 
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8.5s. The compression force reaches its maximum force at 6.1 s and the mold- 
profile switches to force control. The solution to the fields p, Axx, Ayy, Axy is 
plotted as mean values over the mold length. Additionally, the pressure at 
several positions along the flow path is computed in a post-processing step to 
enable comparison with sensors in a press rheometry experiment. The exact 
sensor positions depend on the employed tool and are given for example in 
Section 7.1.4, where the one-dimensional model is used to determine friction 
parameters. 


A simple example without a switch-over to force control is used to illustrate 
the general effect of friction, viscosity and compressibility. The resulting pres- 
sures are plotted in Figure 4.6 and parameters for the frictionless base model 
example are given in Table A.1. For these parameters, the pressure sensors are 
triggered one after another and jump to identical stress levels until the mold 
is filled after triggering the last sensor. The introduction of hydrodynamic fric- 
tion (4=50 N s m3) leads to a pressure difference between sensors, as friction 
introduces a pressure gradient that increases the pressure level from the flow 
front towards the inward direction. An increase of the viscosity (7=2000 Pa s) 
simply shifts the general pressure levels up. The times, at which sensors 
are triggered, remains unchanged in these three configurations, as the bulk 
modulus is chosen high enough to represent quasi-incompressible behavior. 
Lowering the bulk modulus (Kg=1900 Pa) results in a compression of the stack 
and thus a delayed flow front advancement, which can be observed by the 
delayed time at which sensors are triggered. 


Thus, the pressure sensors in the press rheometer deliver an unambiguous 
understanding of the SMC flow. The simple 1D model allows for a quick 
computation of the underlying effects and can be used to fit parameters to 
experimentally recorded results from the press rheometer. 
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Figure 4.5: Exemplary solution of the one-dimensional compressible flow in a press rheometer. 
The maximum compression force is reached after 6.1 s in this example before the filling com- 
pletes at 8.5 s. The different colors in the pressure plot represent the pressures at eleven different 
sensor positions along the flow path. 
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Figure 4.6: Parameter variation in press rheometer model. 
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5.1 Method development 


The fundamental idea of the proposed Direct Bundle Simulation (DBS) ap- 
proach utilizes the observation that fiber bundles often remain in their bun- 
dled configuration during SMC compression molding. Hence, fiber bundles 
can be represented as individual instances in a direct Stokesian dynamics 
simulation with a reasonable number of fiber bundle instances. The num- 
ber of bundles is roughly two or three orders of magnitude smaller than the 
number of fibers and thus computationally accessible. 


5.1.1 Assumptions and scope 


Fiber bundles in SMC can be considered very flexible with a dimensionless 
bending stiffness S* « 1 due to the high aspect ratio and high viscosity of the 
suspending matrix (compare Equation (2.8)). The assembled structure of a 
bundle decreases its bending stiffness even further compared to a homoge- 
neous bundle with the same cross section, because individual fibers may slip 
relative to each other during bending. Therefore, fiber bundles are modeled 
with truss elements, which are one-dimensional elements that only have a 
stretch and compression elasticity, but offer no resistance to bending load. 
Upon bending, individual truss elements remain straight and curvature is 
approximated by angles at the connecting nodes. As discussed previously in 
Chapter 4, the matrix is considered purely viscous and curing is neglected 


79 


5 Direct Bundle Simulation 


during the flow process. The assumptions for the DBS model are summarized 
as 


Fiber bundles remain in a bundled configuration during the molding 
process 


Fiber bundles have negligible bending stiffness 
* The matrix material is homogeneous and purely viscous 


* Curing does not affect the flow process 


5.1.2 Generation of an initial bundle configuration 


First of all, the DBS requires an initial bundle configuration that resembles 
the configuration in the SMC sheets at the beginning of the compression 
molding process. The physical distribution results from bundles falling down 
on the lower foil in a random process (see Section 2.3), thus generating a 
random transversely isotropic in-plane orientation. 


ez E 


e 
Draw Project — ^ Place €x 


Figure 5.1: In a first step, a bundle direction $ is drawn from a uniform isotropic distribution. 


Itis then projected with a second order fiber orientation tensor and normalized to ). Finally, 


the bundle is translated randomly within the initial stack and protruding bundle elements are 
removed. 
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For the simulation, such a distribution is generated by first drawing random 
directions 


py = sin(6) cos($) ey sin(d) sin($) ey + cos(d)e,, (5.1) 


where È € 4(0,2z) and Ó eU(-1,1) with U (a, b) being a uniform probability 
distribution between a and b. This direction is then mapped to a prescribed 
initial fiber orientation state A°, which is often a planar isotropic state in 
SMC, by 

p= pe] 62 
A fiber bundle comprised of several 1D-elements is generated with this di- 


rection and randomly translated within the rectangular domain of the initial 
stack by a shift vector 


AÊ = Afex + A fey + AZe;, (5.3) 


where AX EU (Axmin, AXmax), AV € U (^ ymin; ^ ymax) and AZ € U (AZmin, AZmax)- 
The variables Armin and A£max are the minimum and maximum allowed 


shifts such that at least one node of the bundle resides within the initial stack 
domain. All elements outside the initial stack are deleted. This procedure is 
illustrated in Figure 5.1 and is repeated until the total fiber volume fraction 
equals the prescribed SMC fiber volume fraction. 


The placement does not consider contacts between bundles. These initial 
contacts are resolved during an initial analysis step in the solver Simulia 
Abaqus/Explicit. Overclosures that cannot be resolved are stored as offsets 
and these offsets are resolved during the simulation. 
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5.1.3 Matrix model 


The matrix is described by similar equations as the macroscopic model 


u pikas (5.4) 
p o * div((pv) & v) = div(o) + f” (5.5) 
sn + div((pcy T) -v) = -divid) +0: D. (5.6) 
However, the stress tensor is 
o=-p(p)I +2n(T,y).D' (5.7) 


with an equation of state for the entire composite material and the matrix 
shear viscosity (T, y). The volumetric stress describes the entire SMC com- 
pression behavior including fiber bundles and voids, the deviatoric part 
describes only the matrix response without fiber bundles. A hydrodynamic 
body force field f" has been added to model the influence of fibers on the 
matrix motion and is described in the subsequent section. The body force 
field on the matrix is denoted with a hat symbol (ê) to distinguish it from 
the hydrodynamic force on bundles f}. The evolution equation for the fiber 
orientation state is dropped, because orientation tensors can be computed 
in post-processing from the simulated fiber bundle architecture. 


5.1.4 Bundle-Matrix interaction! 


Bundles are represented as a chain of one-dimensional truss elements (=bun- 
dle segments) that do no transfer bending loads and the stiffness in elonga- 
tional direction is determined by the fiber material. These truss elements do 
not share an interface with the Eulerian domain of the matrix, but float in the 


] Parts of this section are based on [195]. 
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same domain. The fiber-matrix interaction follows the concept of Stokesian 
dynamics (see Section 2.2.6.2) and computes the hydrodynamic force on each 
bundle segment from its velocity relative to the Eulerian environment. The 
resulting hydrodynamic force of each bundle segment f is then distributed 
similarly to the approach by Lindstróm and Uesaka [85] as body force field 
JH on the bulk matrix. 


5.1.4.1 Parameterization of hydrodynamic resistance from micro 
model 


The hydrodynamic resistance of idealized cylindrical bundle segments is 
computed from a numerical study of micro-simulations. The fluid flow in 
a cube is simulated, while a cylinder is placed in its center. The fluid is 
subjected to a unidirectional flow with inlet boundary condition vo; normal 
to the inlet face and a zero-pressure outlet condition, as shown in Figure 5.2. 
All other walls of the cube are allowed to slip, but the normal velocities are 
constrained to zero. A no-slip condition is applied at the cylinder surface. 


The cylinder aspect ratio rp is varied in the range {1,2,3,5,8, 13, 25} and its 
orientation angle relative to the unidirectional flow & is varied in the range 
{0,15°,30°,45°,60°,75°,90°}. The incompressible steady Stokes flow problem 
0=div(-pI+nD') (5.8) 
0=div(v). (5.9) 
is solved for each configuration utilizing the commercial FEM solver Comsol 


Multiphysics. 


To evaluate the resulting forces on the cylinder segment, its surface is param- 
eterized as r = rie, + rte, + rpp, where {e;, e; p} describes the local cylinder 
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Figure 5.2: Two-dimensional section through the cube, in which a Stokes flow is solved to 
compute the drag force on bundle segments. A constant velocity is prescribed at the inlet and a 
zero-pressure condition at the outlet. Other walls of the cube slip, while the cylinder interface is 
subjected to a no-slip condition. The cylinder is parameterized by its orientation angle @ and 
aspect ratio rp. 


coordinate system. With this parameterization, the lateral cylinder surface is 
defined as 


Cu las Torp) ER? In = R,0< 1 < 2,0< ry « ik, (5.10) 
where R is the cylinder radius and / is the length of a bundle segment. The 


surfaces at the cylinder ends are excluded, as these are typically connected 
to the next segment of the bundle and thus do not contribute to the overall 
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resistance. The total hydrodynamic force exerted on the cylinder can be 
determined using an integral over the lateral surface C as 


fil= fonda (5.11) 
C 


with surface normal n. 


A dimensionless measure for the hydrodynamic force in flow direction and 
the transverse direction is obtained with the non-dimensional forms 


1 


= —— — Oy:n dA 5.12 
67t1] R Ve Ee S ; 


ka [024 and k= 
C 


6711] R vao | 
from the vertical and horizontal surface stress components ox and ay. The 
non-dimensionalization is such that the drag factor kq becomes one for a 
suspended sphere and the lift factor kj becomes zero for a suspended sphere. 
Hence, the factors describe the resulting hydrodynamic force components in 
flow direction and perpendicular to the flow direction relative to a suspended 
sphere. 


Drag factor kq Lift factor ky 


~ 
© 
| 


Orientation angle $ in ° Orientation angle $ in ° 


Figure 5.3: Dimensionless drag coefficient k and lift coefficient kj from computation (dots) and 
fit according to Equation (5.13) and Equation (5.14) [195]. 
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The resulting factors from the computation are plotted in Figure 5.3 as dots. 
For small aspect ratios, the factors approach those for a sphere independent 
of the orientation. The drag kg increases with larger aspect ratios and with 
an orientation closer to 90° due to an increased cross section normal to the 
flow direction. The lift component kg peaks at p = 45° and vanishes for a 
fiber aligned with the flow ¢ = 0° or perpendicular to it with @ = 90°. As 
for the drag, the severity of the peak increases with the aspect ratio. The 
computational results are approximated by a fit for the drag factor 


kalrp,®) = 1— a(rp — 1) cos(29) + B(ry — 1) (5.13) 


and the lift factor 
kirp 9) = a(ry — D sin(29) (5.14) 


with a = 0.09 and f = 0.3125. These fits allow efficient evaluation of the 
factors for a given orientation and aspect ratio. They are plotted as solid lines 
in Figure 5.3. 


5.1.4.2 Computation of the relative velocity 


The application of a Stokesian hydrodynamic force requires knowledge of the 
relative velocity between bundle segments and their environment. Figure 
5.4 illustrates how a fiber segment of length / with index j € B at position x; 
and orientation p; floats through the Eulerian domain. The neighborhood of 
such a segment is defined as 


Nj: e£ |0« [Az <n, (5.15) 


where £ describes the set of all indices of Eulerian elements that have an 
element volume fraction e » 0.5. The condition describes a sphere with a 
radius equal to the segment length around the segment center. The neigh- 
borhood relation is dynamic, as fiber bundles move independently from the 
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Eulerian mesh. It has to be updated every time step and the performance of 
the neighborhood search is critical to the overall performance of the DBS. 


Figure 5.4: A fiber bundle with index j and orientation p; is positioned at æ j. The neighborhood 


is called Mj and one exemplary neighborhood element with index i is highlighted at position 
xi. The velocity at this element is v; and its relative position to the bundle segment is termed 
Az; j. [195] 


A naive neighborhood search would include a check of condition (5.15) for 
each bundle and for each Eulerian element. This simple approach is dis- 
played in Algorithm 1 and the time complexity of single run is O(NgNe) 
with the number of bundle segments Ng and the number of filled Eulerian 
elements Ne. 
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Algorithm 1 Pseudo code for naive neighborhood search 


1: for j € 5 do 

2: for i € Edo 
3 Compute Ax; ; 

4 if[Az;;| < 1 then 

5: Compute relative velocity v; — v; for Equation (5.16). 
6: end if 
7 

8: 


end for 
end for 


For the number of elements in consideration here, this algorithm is too slow 
and some additional time to build more sophisticated data structures is neg- 
ligible compared to potential gains during search. Such a more efficient data 
structure is a kd-tree, which is a space-partitioning binary search tree. For a 
1d-tree, i.e. a tree that structures one-dimensional data, the data structure is 
built by sorting the data and taking the median as a node with a left and right 
node containing data that is arranged in the same way. A recursive binary 
bisection allows to find data points very effectively, as each choice for a next 
node rules out approximately half of the nodes at that stage [196]. A kd-tree 
is a generalization, in which a k-dimensional tree is build and searched, as 
indicated in Algorithm 2. The function search is evaluated recursively start- 
ing from a root node and decides which of its two child nodes is closer to 
the target position. If the node farther away is still within potential prox- 
imity of the target value, it is searched, too. Since the function search has 
O(log Ne) complexity due to the effective data structure, the overall perfor- 
mance is O(Nglog(Ne)), which is a massive improvement over the naive 
search method. 


The tree is then used to determine A/; for each segment during a time step and 
compute the relative velocity. The relative velocity for each bundle segment 
is determined as 


Avj= Y, S (vi - v;) (5.16) 
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Algorithm 2 Pseudo code for kd-tree search 


1: Build tree with nodes of € 
2: for j € B do 

3: Nj = SEARCH(root, æ j) 
4 for i € V; do 

5: Compute relative velocity v; — v j for Equation (5.16). 
6: end for 
7 

8 


: end for 


9: function sEARCH (node, x) 
10: if node has no children then 


11: Process terminal node 

12: return results 

13: else 

14: Determine closer node and farther node 
15: SEARCH (closer node, x) 

16: if farther node is closer than / then 

17: SEARCH (farther node, æ) 

18: end if 

19: end if 


20: end function 


with 
2 
_9 [oni | 
wijce ? ? and Wj= } wij. (5.17) 
ieN; 
The weighting factors w;; and normalization with W; represent a Gaussian 
weighting, where the variance is chosen to be //3, such that the sum of 
weights outside V; represents less than 1%. The Gaussian weighting is a 
more accurate scheme than simply taking the nearest Eulerian element for 
the relative velocity and ensures that Av j=vjr-v ls . for — 0. 
J 
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5.1.4.3 Body force field evaluation 


The computed factors ka(rp, od), Kp, œ) and the relative velocity Av j are used 
to compute the force on each bundle segment f" and the corresponding 
body force field f". The hydrodynamic force on each element is 


fË = GART} | kat. G)Avj+ hp.) |Avj] aj] 6.18 
with the angle between the relative velocity and the fiber orientation 
Avj: Pj 
Qj = arccos Pi (5.19) 
[sv 
and with the direction perpendicular to v; computed by 
qj = -sgn(p; :Av;) (er [av]} lav] (5.20) 


in the plane that is formed by the two unit vectors p; and [av il (details may 
be found in Appendix A.1). 


Finally, the same Gaussian weights from the determination of the relative 
velocity are used to distribute the body force contribution of each bundle. 
Hence, the contribution of bundle j to the body force field in Eulerian element 
iis 

apo 1 Wij 


Hh (5.21) 


with the volume of the Eulerian element V;. These contributions are summed 
up for each Eulerian element. 
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5.1.4.4 Implementation 


The model is implemented with a combination of several subroutines in the 
commercial FEM program Simulia Abaqus/ Explicit. 


A VEXTERNALDB subroutine is utilized to control the entire analysis workflow. 
The subroutine parses the input deck before the start of the analysis to extract 
parameters such as index ranges of bundle elements and Eulerian elements. 
During the analysis it builds kd-trees at each time increment for all elements 
that have material volume fraction e > 0.5. The implementation for building 
and searching the kd-trees itself is taken from Kennel [196], who built one of 
the most efficient kd-trees in Fortran. A tree is build at the beginning of each 
time step using new memory such that the previous tree is still in memory 
for delayed parallel threads working on this tree. The memory is freed such 
that always two alternating trees, an old one and a new one, exist at all times 
(see Figure 5.5). Additionally, the VEXTERNALDB subroutine computes an 
approximate stable time increment for the bundles 


R 
Atg = 0.25 


(5.22) 


Amax 


based on the idea that a bundle should travel no more than a fraction of a 
bundle radius during a step with the maximum acceleration of all bundles 


Amax: 


A VUFIELD subroutine is used to extract positions and velocities at nodes 
at each time step. The values are saved to field variables to make them 
accessible to a VUSDFLD subroutine, which interpolates position and velocity 
at integration points of Eulerian elements as well as truss elements. Each 
Eulerian element (EC3D8R or EC3D8RT) and each truss element (T3D2) has 
only one unique integration point. The variables e, T,V,n,x,v, fH fE at the 
integration points are saved to global arrays with a length equal to the sum 
of truss and Eulerian elements for further processing in other subroutines. 
The global arrays are accessed with the index of elements and due to the 
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Time 
Å] 


Memory 
Tree 1 


First tree remains in memory Overwrite with third tree 


Memory 
Tree 2 


Build first tree Build second tree Second tree remains in memory 


Figure 5.5: The first tree is built using allocated memory for the maximum tree size at the 
beginning of the first time increment. This tree should not be overwritten at the beginning of the 
next time step, as some threads may still operate on this data structure. Hence, the second tree is 
built in different pre-allocated memory. Then, the third tree overwrites the memory of the first 
tree to save time on allocation/deallocation and the second tree remains in the second allocated 
memory space. This procedure with alternating trees omits any waiting for threads, which is 
potentially disastrous for threaded parallel performance. 


parallel domain decomposition, the access happens in parallel while still 
being thread safe. 


Finally, a VDLOAD subroutine evaluates equations (5.16), (5.18) and (5.21). It 
performs parallel tree searches in the current kd-tree, which do not alter the 
tree, but only obtain the neighborhood relation. 


5.1.5 Bundle-Bundle interaction? 


The objective of this section is a model for short-range interaction of bundle- 
bundle contact points. The model needs to determine the mechanical contact 
in normal direction and tangential friction as well as lubrication in normal 
an tangential direction. 


A simple mechanical contact can be modeled quite easily as repulsive force 


0 >0 
fas = si (5.23) 
Kchap g <0 


2 Parts of this section are based on [32]. 
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with a contact stiffness K. and the gap between bundle surfaces g. The force 
is absent, if fibers are separated by a gap, and penalizes overclosure linearly. 
The corresponding tangential friction force is 


; 0 g>0 
fap [Avas] g<0 


B7 -u| D (5.24) 


with a Coulomb friction coefficient u. The friction acts opposed to the direc- 
tion of the relative tangential velocity Ava. 


The Newtonian form of the tangential lubrication given in (2.73) is 


d? Av 
ny—— arg 
| Pa XP 
fi = ß (5.25) 
contact area effective shear rate 
0 gs0 


and other models [30, 56,57] may be expressed with it using a proper choice 
of the effective sheared gap G. This work uses the effective sheared gap as a 
single adjustable parameter, because it can be interpreted as the mapping 
from a physical distance between bundles g to an equivalent sheared effective 
gap G accounting for all geometrical details of the contact area. This relation 
is obtained in the following two sections using analytical considerations and 
a parametric simulation study. 


5.1.5.1 Analytical considerations for the effective sheared gap 


The effective sheared gap G is equal to the separation distance g in case of 
two flat fiber bundles with parallel faces. In this case, the interaction force 
becomes infinite, if the separation distance between two such bundles ap- 
proaches zero. This section aims at finding an expression of G for bundles 
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with elliptical cross sections and investigates the behavior when the separa- 
tion distance becomes close to zero. 


Bundle a dal sin(y) 


Figure 5.6: Schematic illustration of the contact area between two fiber bundles. The first bundle 
a is aligned with the x-direction and forms the angle y between this bundle and a second bundle 


B. 32). 


Two identical fiber bundles with constant elliptical cross sections are consid- 
ered. The bundle cross section is described by a minor axis diameter d, and 
a major axis diameter d4. Both bundles are positioned with their minor axes 
aligned in z-direction, as depicted in Figure 5.6. Bundle a is aligned with the 
x-direction and a second bundle f crosses it such that an angle y is formed 
between them in the xy-plane. Bundle f is placed such that a separation 
distance g remains at the closest contact point between both bundles. The 
coordinates (x, y, z), the separation distance g and the effective sheared gap 
G are non-dimensionalized with the bundle dimensions as 


* 


Haw] ll, pf pof, 0-8 — ga 
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to obtain generalized results. The dimensionless z*-component of bundle 
surface a may be expressed as 


1 
a=- Ayri- (5.27) 


with y* € [0,1]. The dimensionless z*-component of bundle surface ß is 
consequently 


1 1 


£57 ~ "E -(x*- y* cosy)? -g' (5.28) 


with x* € -i + y* cos(y), i 4 y* cos(y)]. The upper and lower bound de- 
pend on the y-position as this bundle is placed with an angle y relative to 
bundle a. The separation distance is added to the parameterization of this 
bundle surface to place the point closest to the surface of bundle a at distance 
g relative to bundle a. 


A strong simplification is introduced by assuming that the velocity profile 
between both surfaces is approximately linear. Hence, the velocity gradient 
at Pus infinitesimal surface element of the interaction zone is given as y — 


Ê with the velocity difference of the bundles Av, p. Assuming Newtonian 


ZB Za 
viscosity 7, an integration over the sheared domain yields 
d.d? I y" cos(y) 1 
Fi, topi a [ ts —LLdx'dy. — (629 
|sin(y)| —}+y* cos(y) 25 -2g 


1/G* 


Inserting equations (5.27) and (5.28) and using the substitution s* = x* — 
y* cosy + 3, this can be simplified to 


1 
= ds* dy", (5.30) 
I 1- Vs*(1-s*)- yy” (1 -y*) g* 


where it becomes apparent that the result is independent of y. 
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Solving the inner integral leads to 


LET | Akh(gy) - 


- (3 eme ani] (5.31) 
y 4kilg*, y*)*- 


with kı =1+g*-,/y*(1- y*). After a change of integration variables, with an- 


other substitution kz = yag”? —8g*u* +4u** + 8g* — 8u* +3 and the first- 
order approximation arctan x = 7x for small values of g*, this is further 


simplified to 


(5.32) 


i af u[(2+g* - u*)wke(g*,u*)—- 4g" —4u" n] | 
G* 0 ko(g*,u*)V—4u** +1 


Unfortunately, the result of this integration can only be expressed in terms 
of incomplete elliptic integrals. However, a series expansion leads to an 
analytical expression that is exact if g* approaches zero. Employing only the 
first series term, the effective sheared gap for two bundles can be expressed 
as 


G* = 2 -O(g*) (5.33) 


n(-16arctanh 2 + 16/3 - 48 + 241n2 - 8lng* +7) 


The dimensionless effective sheared gap depends only on the dimensionless 
separation distance g* and is independent of the contact angle y and the 
elliptical cross section. The effective sheared gap approaches zero, if the sep- 
aration distance approaches zero. This means that interaction forces would 
become infinite for infinitesimal close contact, just as in the case of flat bun- 
dles. However, the assumption of linear velocity profiles is strong and should 
be compared to full solutions of Stokes' equations on the sheared domain. 
Therefore, the next section presents a parametric study of the sheared gap 
that solves the creeping flow between two bundles numerically. 
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5.1.5.2 Parametric study of the effective sheared gap with a micro 
model 


The sheared fluid between two fiber bundles is simulated to determine the 
effective sheared gap t between bundles. A parametric model of the sheared 
domain is built using the following parameters for a parametric sweep: 


Aspect ratio of the cross section The aspect ratio of the cross section is de- 
fined as r. = d4/dy. The minimum value in this application is one, 
which corresponds to a circular cross section. The maximum value 
investigated here is re = 8. 


Contact angle The contact angle between two bundles is described by y. 
This value ranges from 15° to 90° in 15° steps, where 90° corresponds 
to perpendicular bundle orientations. The angle 0° would lead to a 
singularity and is excluded. 


Separation distance The separation distance g describes the distance be- 
tween the surfaces of both bundle geometries at the closest point. The 
gap cannot be reduced to zero, because this would form a singularity 
due to the infinite shear rate at the contact point and would result in 
poor element quality in the sheared domain. 


Load angle The direction of relative bundle motion between bundles is de- 
noted as load direction 6. Due to symmetry, 6 is varied between -90° 
and 90° in 15° steps. 


A small buffer zone around the domain is added to the sheared domain to 
prevent poor element quality at the sharp boundaries. An illustration of the 
parametric model for one exemplary configuration is given in Figure 5.7. 


Stoke's equations for creeping flow are solved on the domain highlighted 
green in Figure 5.7. The momentum balance is 


0 = div [- pI+n (grado) + grado)” )] (5.34) 
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Figure 5.7: Two fiber bundles with a sheared domain (green) between them. The elliptical cross 
section is characterized by minor and major diameter dp and da, respectively. The separation 
distance between bundles is denoted as g and the upper bundle is placed with an angle y 
towards the lower bundle. The load direction rn (=direction of relative motion) is characterized 
by an angle 6 ranging from -90° to 90° due to symmetry. [32] 


with pressure p, Newtonian viscosity 7 and velocity v. The fluid is assumed 
incompressible, hence the mass balance gives 


div(v) = 0. (5.35) 


The boundary conditions are illustrated in Figure 5.8. The wall in contact 
with bundle a is fixed with zero velocity. The wall in contact with bundle 
p is subjected to a tangential velocity with magnitude vo and direction 6 
within the shear plane. The remaining walls are subjected to a zero-pressure 
boundary condition and the inflow or outflow of fluid is limited to a direction 
m that corresponds to the load direction. 


Figure 5.9 visualizes some computed effective sheared gaps G* for different 
aspect ratios of the cross section in each sub-figure. Each set of lines indicated 
by the same plot style is composed from simulation results for contact angles 
y € [15°,30°, 45°, 60°, 75°, 90°]. Conversely to the analytical solution, these 
results vary with load direction and small peaks occur if the load direction 
matches the orientation angle. 
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Fixed walls Moving walls Pressure and direction 
Se uL « — 
vo cosó cosó 
v=! 0 v-| wsind m-| sind 
0 0 


Figure 5.8: Boundary conditions for a sheared volume with rc = 2 and $ = 45°. The blue areas 
indicate the surfaces at which boundary conditions are applied. The vectors v denote the velocity 
boundary condition at the corresponding surfaces. The vector rn denotes the constrained 
direction of velocity at the Neumann boundaries and depends on the load angle 6. [32] 


Figure 5.10 shows the mean value of the dimensionless effective sheared gap 
G* plotted against the dimensionless contact gap g* for all contact angles 
and load directions. The colored lines represent the simulation results for 
different aspect ratios of the cross section, which correspond to the colors in 
Figure 5.9. The error bars indicate the standard deviation within each set due 
to varied contact angles and load angles. The analytical result of Equation 
(5.33) is plotted as a solid black line. 


The mean values of the simulation results fall on one curve in Figure 5.10 
regardless of the cross section aspect ratio rc, even though the models vary 
significantly in size and shape. The simplified analytical solution seems to 
overestimate the effective sheared gap. This may be attributed to the pressure 
distribution and resulting non-linear velocity profiles across the gap, while 
the velocity profile is assumed linear in the analytical solution. 


Literature values for the effective gap G* are given for example by Le Corre 
et al. [8] and Guiraud et al. [104]. Le Corre et al. estimate a value G = 
2.0 x 10? mm from compression experiments for bundles with minor di- 
ameter dp = 0.06 mm, which is equivalent to G* = 0.033. Guiraud et al. report 
G = 1.5 x 10? mm from pull-out experiments and dp = 0.2mm, which is 
equivalent to G* = 0.075. Both literature results are added to Figure 5.10 
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Figure 5.9: Dependence of G* on load angle for exemplary simulation results with varying 
contact angles y (line style) and varying aspect ratios (color) at different separation distances g*. 
Not all simulations are included to enhance visibility. Mean values of all load angles and contact 
angles are used to generate the diagram in Figure 5.10. [32] 


at g* = 0 for reference. The magnitude of these reported values is in agree- 
ment with the results of the analytical considerations and parametric study 
presented here. 


The analytical considerations showed that the interaction force between two 
elliptical bundles becomes infinite for infinitesimal small separation distance. 
However, this singularity cannot be introduced to numerical simulations. 
Therefore, an offset for the effective sheared gap Gj ~ 0.04 is proposed and 
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Figure 5.10: Comparison of several approaches to determine the dimensionless effective sheared 
gap G* asa function ofthe dimensionless separation distance g*. [32] 


used in a simple fit of the results shown in Figure 5.10. One possible fit for 
the simulation results is given as 


G* = 0.027arctan(65.5g") + g* + Gy. (5.36) 
This approach fulfills the condition 
G* |g*—0 = Go (5.37) 
to achieve a finite separation for bundle contact and the condition 
G*|g* oo = g” (5.38) 
for large separations where the surface curvature becomes irrelevant com- 


pared to the separation. 
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Figure 5.11: The variables passed into the subroutine (here g) assume a cylindrical shape (dashed 
lines), which overlap in this case. However, bundles of thickness dp may not overlap and a gap 
remains. 


5.1.5.3 Implementation 


Bundle-bundle interactions are implemented with a VUINTERACTION sub- 
routine in Simulia Abaqus/Explicit. By default, the contact variables of any 
truss or beam contact imply a cylindrical shape, but one may modify these 
parameters in the subroutine to account for a different shape. In this case, 
the gap is defined as 

g-d4-d, (5.39) 


where d = \/dady is the equivalent cylinder diameter, is the gap between 
cylinders as provided by Simulia Abaqus/Explicit' contact solver and dp is 
the bundle thickness. This definition is a simple approximation to correct 
the contact thickness from a cylindrical surface to an elliptical shape, as 
illustrated in Figure 5.11. Strictly, this correction applies only for bundles that 
are perfectly aligned with their minor axis and neglects contact of bundle 
edges along the major axis. However, a contact of bundle edges is rare in SMC 
compared to the contact of the relatively flat bundle surfaces. Further, the 
rheology is likely dominated by the large sheared zones between bundles, 
which share a contact area in the plane of their major diameters. 
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In this implementation, the normal contact stress is defined as 


—Keg+KaAımn g<0 


Sn = di Av (5.40) 
Tl siny| GA. g>0, 


where Ke is the contact stiffness (typically Young's modulus), Kg is a contact 
damping parameter (typically 7), Ac is the contact area and Avy is the rela- 
tive velocity of the contact partners in the normal direction. The tangential 
contact stresses are defined as 


-u| fe; | [Av] g <0 
8 uw (5.41) 
~Iisiny] GAc g>0 


where Av; is the relative velocity of the contact partners in tangential direc- 
tion. The dependence of G on the contact gap g is modeled with Equation 
(5.36). 


5.1.6 Post-processing of bundles 


The DBS results in a bundle configuration represented by a set of nodes and 
line elements for each simulation time step. Post-processing is required to 
derive macroscopic properties of the result and compare it to macroscopic 
models and experimental results. 


5.1.6.1 Bundle curvature 
The bundle curvature is a measure for buckling of fiber bundles and is typ- 


ically high in areas that experience compression during the process. The 
curvature is approximated as 


2 1 
Ck = 1 tan 2 arccos (pi 222 (5.42) 
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at each node k connecting two neighboring bundle segments of length / 
with unit directions p; and pj. This evaluation assigns each nodal position 
an approximate curvature value, which can than be compared to curvature 
results from computer tomography. 


5.1.6.2 Contact number 


To compute the total number of contacts between bundles, a distance map 
between bundle element centers is computed. This distance map describes 
for each bundle element center the distance to each other bundle element 
center. Potential contacts occur, where the distance between two centers j 
and j' is smaller than the segment length /. For these possible contact points, 
the normal to both directions is 


nz [p;x»;] (5.43) 


and the relative position can be determined by solving the linear system of 
equations 
[p;x Nx -py.s| [ui] 
wj-aj=|pjy ny -pjy| uo (5.44) 


Pjz Nz Piz u3 


for u to ua. A contact occurs, if| u| < 1/2 and | us| < 1/2 (that is, the intersec- 
tion is within the range of both bundle lengths) as well as | up| < d (that is, the 
intersection is close enough to be considered a contact of overlapping cylin- 
ders). The evaluation of contacts during the simulation allows to estimate the 
influence of bundle-bundle interactions on the resulting compression force. 
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5.1.6.3 Mapping of discrete orientation tensors and fiber volume 
fraction 


If the fiber architecture has to be compared to a macroscopic model or if the 
results of the process simulation are required as input for further structural 
simulation, a mapping of line elements to an arbitrary mesh is necessary. 
Hence, for each cell i, the length of each potential bundle segment j in this 
cell Al; ; is computed. If bundles cut a face of the cell, the length within the 
cell up to the intersection point is computed, as indicated in Figure 5.12. 


<< ——— HÀ — IT 
Al;; =0 
(a) Both nodes within cell (b) Both nodes outside cell (c) Partial intersection 


Figure 5.12: Case distinction for computation of the length AJ; j in a cell. 


The local fiber volume fraction is then approximated as 
= (5.45) 


which assumes that the remaining cell volume is filled completely by the 
composite. The local discrete second order fiber orientation tensor is 


1 
Se 0 A bs (5.46) 
M jeg; Alij 2x. d 1 


i 


105 


5 Direct Bundle Simulation 


where 5; c B is the set of bundles with a length AJ; ; > 0 in cell i. The proce- 
dure is implemented as a ParaView filter and allows generic mapping between 
CAE files in the VTK file format. ? 


5.2 Verification 


5.2.1 Orientation evolution of a single bundle^ 


The motion of a single bundle in shear flow must be faithfully reproduced to 
qualify the method for computation of fiber reorientation. Therefore, a fiber 
bundle with a length of 25 mm and an aspect ratio of ry — 25 is subjected to 
a shear rate y = 10s! in this verification. The domain for this simulation 
is illustrated in Figure 5.13. The bundle is placed at the center, discretized 
with ten segments and positioned vertically, such that the initial orientation 
is 0 — 0. 


80 mm 


-1 


40 mm 
a 
x 


Velocity in mm s 


Figure 5.13: The contour plot shows a fiber bundle discretized with ten segments in a shear flow. 
The color codes indicate the velocity in x-direction. The fluctuations at both fiber bundle ends 
show how the two-way coupling influences the macroscopic velocity field [195]. 


3 The filter is available at https : //github.com/nilsmeyerkit/paraview map.lines 
4 Parts of this section are based on [195]. 
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Figure 5.13 shows bundle position and velocity in x-direction shortly after 
starting the simulation. The contour plot of the horizontal velocity compo- 
nent depicted in Figure 5.13 indicates the two-way coupled nature of the 
presented approach, particularly visible at the two ends of the fiber bundle. 
Although the bundle is flexible, it behaves like a rigid body until alignment 
with the flow due to the positive normal stress in the direction of the bundle 
axis. 


Figure 5.14 compares the orientation evolution of the DBS with ten truss 
elements and two truss elements to the solution of Equation (2.22) with 
aspect ratio an equivalent aspect ratio computed from the formula by Cox [13]. 
The simulation is in good agreement with the reference solution for both 
discretizations. 
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Figure 5.14: Comparison of bundle orientation angle computed from DBS and Jeffery's Equation 
with an bundle aspect ratio Tp = 25 [195]. 


There is a small difference between simulation and analytical solution at 
the almost horizontal state in Figure 5.14. At this point, torque induced by 
friction at the lateral surface dominates bundle motion. In SMC, bundles 
are heavily confined by other bundles and the mold. It is assumed that the 
torque that spins a free bundle in a dilute situation is small compared to the 
confinement effects and therefore is neglected here. 
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Additionally, a bundle with a bundle aspect ratio rp = 25 is placed 90° to the 
flow under the same conditions as in the parameter identification (see Figure 
5.2 in Section 5.1.4.1) and meshed with one and ten segments. The resulting 
drag force normalized with 67:17] R vos is 9.31 and 9.55, respectively. This is 
close to each other, but slightly smaller than the drag coefficient shown in 
Figure 5.3, because the averaged velocity around the bundle is smaller than 
the nominal velocity far away. Nevertheless, the orientation result and the 
total drag indicate that bundle motion is generally only slightly affected by 
discretization. The choice of bundle discretization may be considered the 
choice of which length scale of velocity disturbance is computed by the 
Eulerian mesh and which length scale is included in the drag coefficients. 


5.2.2 Anisotropic flow front progression 


The presence of fiber bundles leads to anisotropic flow and the extreme 
case of unidirectional orientation is considered in this verification. There- 
fore, an initial rectangular stack of the dimensions 25mm x 25mm x 15mm 
and fiber bundle length 10mm is compressed with 1mms~! between two 
plates and simulation results of the macroscopic reference model, the DBS 
and an isotropic macroscopic model are qualitatively compared. Additional 
simulation parameters can be found in Table A.2. 


Tables 5.1 and 5.2 show the flow front progression of the suspension for 
0.25% fiber volume fraction and 2.5% fiber volume fraction, respectively. 
The left columns show isotropic results, where the stack is evenly stretched 
in horizontal and vertical directions. There is a slight tendency towards a 
more circular shape, as the edges of the initial stack are rounded during 
compression. 


The anisotropic models consider unidirectional horizontal fiber bundle align- 
ment. For 0.25% fiber volume fraction, there is a slight preference for flow 
perpendicular to fiber bundles, as the elongation in fiber bundle direction 
is prohibited by quasi-inextensible fiber bundles. The DBS shows rounder 
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Table 5.1: Semi-dilute anisotropic flow with unidirectional fiber orientation and 0.25% fiber 
volume fraction. 


Time Macro Isotropic Macro Anisotropic DBS 


2.58 


5.0s 


corners than the mascroscopic model, which is probably due to the inhomo- 
geneous dispersion of bundles in the domain (there are just 125 bundles in 
the model) and specifically due to the absence of long bundles in the cor- 
ners, as they are cut outside the stack. In contrast, the macroscopic model 
assumes equal fiber aspect ratios at every material point without considering 
the stacks boundaries. The flow front shape is similar for both models and 
the anisotropy is described qualitatively correct. For a higher fiber volume 
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Table 5.2: Semi-dilute anisotropic flow with unidirectional fiber orientation and 2.5% fiber 
volume fraction. 


Time Macro Isotropic Macro Anisotropic DBS 


"n i j J ig) 


2.5s 


5.0s 


fraction of 2.5% the anisotropy becomes more severe and there is much more 
elongation perpendicular to fiber bundles than in fiber bundle direction. The 
round corners of the DBS model likely arise from a lack of long fibers in the 
corners. The overall comparison shows that DBS can describe the anisotropic 
motion of the fluid qualitatively correct. 
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5.2.3 Frictionless linear compression 


To verify quantitatively correct compression forces due to anisotropic viscos- 
ity, a quasi-incompressible press rheometer configuration (see Figure 4.4) 
without friction is considered in this section. Therefore, an initial rectangular 
stack of the dimensions Xo = 50mm, B = 50mm, ho = 5mm is compressed 
with a constant compression speed hy = 1.66 mms! between two rigid plates. 
The mold length is Xmax = 75 mm, which results in an initial mold coverage 
of 66% and a volumetric fill time of 1 s. Snapshots of the mold filling process 
of such a configuration with 25% fiber volume fraction using DBS are shown 
in Figure 5.15. 


t=0s 


t=0.4 s 


t=0.8 s 


Figure 5.15: Snapshots of the compression molding process with 25% fiber bundles. The top 
mold is not displayed for visualization purposes. 
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First, the resulting compression force of the pure matrix without any fibers is 
computed by the analytical Equation (4.49), the 1D PDE model given in Sec- 
tion 4.2, the 3D VUMAT model 4.1 and a DBS without any bundles. The results 
are depicted in green colors in Figure 5.16 and show congruent compression 
forces. The compression force increases non-linearly, as the magnitude of the 
compression strain rate D = ho/h increases non-linearly, since h is reduced 
with a constant rate. The results at the very beginning are left out of this 
visualization for all numerical models, because the material is initially at rest 
and has to accelerate after establishing a contact with the mold. 


To verify the viscosity increase due to added bundles, 300 bundles with a nom- 
inal length of L= 25mm are added employing the procedure described in 
Section 5.1.2. After the removal of protruding bundle elements, this equates 
to approximately 1% fiber volume fraction in the stack. The DBS model 
is then solved with disabled bundle-bundle interactions and with enabled 
bundle-bundle interactions assuming a bundle width dj = 0.5mm. As indi- 
cated by orange bullets in Figure 5.16, the introduction of 1% fiber bundles 
approximately triples the reaction force compared to pure matrix. The in- 
corporation of short-range bundle-bundle interactions hardly changes the 
result, as interaction between bundles is relatively rare. On average, the 1610 
bundle segments experience 131 contact pairs resulting in only 0.16 contacts 
per bundle segment, or in other words, every sixth fiber bundle is in potential 
contact with another bundle at all. 


For reference, the analytical solution, 1D model and 3D VUMAT model are 
solved for 1% fiber volume fraction as well. These models require an estima- 
tion of a fiber aspect ratio that is comparable to the generated architecture 
with a length distribution due to the cutting process. The average fiber bundle 
length is computed as 


jn-t pp dxd 
E= Í r L(x, y)dx (5.47) 
XoBo Jo Jo 2 » 
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10° 
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Compression force Fc in N 


d — Analytical (f = 0%, kp =0) +-+- 1D model (f = 096, kp = 0) 
=== VUMAT (f = 096, kp = 0) @ DBS(M) 

102 -| —— Analytical (f = 1%, kp 20) "+ 1D model(f = 1%, kp = 0) 

| === VUMAT (f=1%, kp =0) —— VUMAT (f = 196, kp — 5) 

| € DBS (M+BM) E DBS (M+BM+BB, 0.5mm) 

J — Analytical (f = 25%, kp = 0) === 1D model (f = 25%, kp = 0) 

| --- VUMAT (f = 2596, kp = 0) @ DBS (M+BM) 

- VUMAT (f = 2596, kp =5) DBS (M + BM + BB, 0.25 mm) 

DBS (M + BM + BB, 0.5 mm) DBS (M + BM + BB, 1.5 mm) 

10! j j j | j j j j | 

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 


Compression time t in s 


Figure 5.16: Comparison of analytical, 1D macro, 3D macro and DBS models to verify the 
implementations and models. The abbreviation M, BM, BB stand for a matrix component, two- 
way-coupled bundle-matrix interaction and bundle-bundle interactions, respectively. 


and results to L = 7/12L in this case, because the initial stack is small in 
comparison to the bundle length. The aspect ratio definition is not unam- 
biguous for fiber bundles and candidates are the ratio between some average 
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fiber length and fiber diameter, equivalent bundle diameter, major bundle 
diameter. For simplicity, the aspect ratio between average fiber length and 
equivalent bundle diameter is chosen here, which is rp = 72.9. The analytical 
solution depicted as solid orange line in Figure 5.16 uses the constant values 
Axx = 0.5, Ax = 0.375 for the fiber orientation state, which represents an 
ideal planar isotropic distribution. The solution of the 1D model is depicted 
as dotted orange line and includes the computation of fiber orientation evo- 
lution. It starts at the same compression force level as the analytical equation, 
but the force increase exceeds it, as fibers align with the flow. The 3D VU- 
MAT solution (orange dashed) matches the 1D solution precisely for disabled 
short-range bundle-bundle interactions. With bundle-bundle interactions 
(kp = 5, solid light orange line), the compression force is increased by a rather 
small amount, which is in line with the expected behavior at this low fiber vol- 
ume fraction and with the increase introduced by bundle-bundle interactions 
in DBS. 


At 25% fiber volume fraction, the compression force level of DBS without 
bundle-bundle interactions is significantly increased compared to the previ- 
ous results due to the increased in-plane elongational viscosity. Again, the 1D 
solution starts at the same compression force level as the analytical solution, 
but deviates with increasing time due to fiber alignment (violet solid line and 
dotted line). The 3D VUMAT solution shows a slightly higher compression 
force response, as this simulation requires usage of a shear number N, = 10 
to stabilize the numerical solution procedure. Very high ratios of elonga- 
tional viscosity to shear viscosity tend to cause transverse shear bands that 
destabilize the solution process. 


In the DBS model with 25% fiber volume fraction, 87120 contacts pairs form 
between the 37577 involved bundle segments on average, which means that 
each segment is in contact with 4.63 other bundles. Therefore, a signifi- 
cant contribution of short-range bundle-bundle interactions is expected. 
Hence, bundle-bundle interactions are introduced with 0.25 mm, 0.5 mm 
and 1.5mm bundle width in the DBS model and indicated with light violet 
stars, squares and triangles in Figure 5.16, respectively. The cross section area 
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of bundles is not modified, such that wider bundles are proportionally flatter. 
It becomes apparent, that at such a high fiber volume fraction, short-range 
hydrodynamic bundle-bundle interactions indeed lead to a significant in- 
crease of the compression force. The compression force increases with flatter 
and wider bundles. For a width of 0.5 mm, the hydrodynamic interaction 
factor kp = 5 yields similar results to the DBS. 


0.7 
E —— One-way coupled BM Axx 
B - - - One-way coupled BM Ayy 
E —— Two-way coupled BM Axx 
8 06 | ... Two-way coupled BM Ayy 
" 
2 —— Two-way coupled BM+BB Axx 
E | ==- Two-way coupled BM + BB Ayy 
8 054 II 
3 ELLA ce 
& ""- 
= "" "mene 
o 

0.4 

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 


Compression time t in s 


Figure 5.17: Evolution of discrete orientation tensor components during the linear compression 
example. The abbreviations BM and BB stand for bundle-matrix interaction and bundle-bundle 
interaction, respectively. 


While the compression force can be altered significantly by the presence 
of bundle-bundle interactions, the orientation evolution experiences only 
little change. This is demonstrated in Figure 5.17 for the DBS with 25% fiber 
volume fraction, where one-way coupling refers to an application of Equation 
(5.18), but not Equation (5.21), which results in the same compression force as 
the pure matrix behavior. The globally evaluated discrete second order fiber 
orientation tensor components are hardly influenced by two-way coupling 
and bundle-bundle interactions. 


This section verified, that the two-way coupled bundle-matrix interaction 
increases the compression force significantly. Without bundle interactions, 
the analytical solution, 1D model, 3D VUMAT model and DBS predict similar 
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orders of compression force magnitudes. The effect of bundle-bundle inter- 
action further increases the compression force. Bundle-bundle interaction 
introduces locally large interaction forces between bundles, which affect the 
time step of the explicit DBS negatively. While bundle-bundle interactions 
have a large effect on compression force, they have a weak effect on fiber 
re-orientation, as many bundle-bundle interactions likely cancel each other 
w.r.t the overall motion of a bundle. If the fiber architecture is of primary 
interest in a simulation, bundle-bundle interactions could be omitted for the 
benefit of better computational times. 


However, the choice of rp and kp is ambiguous. The choice of rp = 72.9 is 
reasonable, but it cannot be considered the definitive result, because the ref- 
erence model assumes rods with constant length, while it is applied here for 
bundles with a length distribution. Additionally, the application of Shaqfeh 
and Fredrickson’s model at 25% is certainly beyond the semi-concentrated 
regime and the aspect ratio is not uniform for all fiber bundles. A validation 
with experimental results of a press rheometer is necessary to validate the 
results and evaluate the importance of contributions from two-way bundle- 
matrix interaction, bundle-bundle interaction and friction at the mold sur- 
face. This will be addressed in Section 7.1.4 and Section 8.3. 
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6.1 Method development 


The patches are local continuous carbon fiber reinforcements that are pre- 
formed to a desired shape and subsequently co-molded with SMC. They can 
be placed in critical areas of a part to provide a tailored improvement of 
mechanical performance while keeping the amount of costly unidirectional 
carbon fibers at a minimum. This chapter addresses the formulation of a 
hyperviscoelastic model including a simple damage model. 


6.1.1 Assumptions and scope 


Typical shell formulations of UD fabric forming assume decoupled mem- 
brane and bending behavior of the fabric, to account for the relative slip 
between fibers during bending [171]. However, the investigated material is 
much stiffer than fabrics in wet compression molding or thermoforming at 
process conditions and does not show a clear separation between membrane 
and bending properties. This assumption is verified in the characterization 
experiments in Chapter 7. Given the high stiffness and the initial temperature 
of the material in the process, the preform’s ability to deform further during 
the molding process is limited. Patches rather disintegrate upon large defor- 
mations instead of being draped into a complex shape during co-molding. 
Consequently, the main objective of the co-molding simulation in this context 
is a prediction of the damage initialization in patches under certain process 
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conditions. Hence, a simple Hashin criterion is implemented to describe 
patch damage by matrix tension during co-molding. 


The patch material considered here is made from a stitch-bonded unidirec- 
tional carbon fiber fabric. It is stitched with glass fiber yarn as sketched in 
Figure 6.1. As a simplification of the stitching pattern, it is assumed that the 
stitching can be represented by a single direction that is initially transverse to 
the carbon fiber. 


SE EEEEEEINIKKKKKK, 
EERS 
HERSES 
mpm AAA AAA 


(a) Front (b) Back 
Figure 6.1: Sketch of the stitching pattern of the fabric used for fabrication of the unidirectional 
reinforcement patches (ZOLTEK PX35 UD300). 


Considering the observed material behavior, the assumptions are summa- 
rized as follows: 


* Membrane and bending do not require a decoupled model 

* The material can be considered hyperviscoelastic 

* AHashin damage model is sufficient to describe the observed damage 
* The patches are thin and can be modeled with shell elements 


* The stitching direction is initially perpendicular to the carbon fiber 
direction 
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6.1.2 Hyperviscoelastic model 


The material response (membrane, bending, shear) is modeled based on an 
Ideal Fiber Reinforced Material (IFRM). The IFRM model assumes inextensi- 
ble fibers and incompressibility, which form a reaction term. The stress can 
be expressed as 

o =—pI+ Sp(pp® pp) +7 (6.1) 

[acne 

with the hydrostatic pressure p, a fiber stress Sr in fiber direction pr and an 
extra-stress T, which is modeled via a constitutive equation [197, 198]. 


Following Dórr [199], a hyperviscoelastic model based on the IFRM is formu- 
lated in the initial configuration as 


S = Se + Syt Sr (6.2) 


with S denoting the second Piola-Kirchoff stress tensor. It comprises an 
elastic (e), viscous (v) and reaction (r) part. The isotropic elastic part is a 
hyperelastic St. Venant-Kirchoff material with 
EuG " 
Se =| — HP, +2GmP2|: E, (6.3) 
3Gm - Em 

where Ey and Gy are the isotropic matrix Young's modulus and shear modu- 
lus, respectively. The notation $ denotes that these properties depend on the 
matrix damage variable. The viscous part is 


S,-Fi.|2. F'-E-F! |.-F", (6.4) 
0 


current configuration 


—Á————— 
initial configuration 
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where the rate of deformation E is transformed to the current frame, then 
multiplied by twice the patch viscosity 7]; and then transformed back to the 
initial configuration. The reaction term is expressed as the regularization 


Sr = Epep (pr e pr) with er=E: (pr ® pr) (6.5) 


representing the extension in fiber direction. Finally, the second Piola- 
Kirchoff stress is transformed to the current configuration by Equation (2.83). 


6.1.3 A simple damage model 


Large deformations of the unidirectional patches are often accompanied by 
damage in this application. Therefore, the prediction of damage initialization 
during the molding process is critical to ensure that patches are designed 
such that they stay intact. 


Hashin developed a model specifically for unidirectional fiber reinforced 
composites and distinguished the four failure modes matrix tension, matrix 
compression, fiber tension and fiber compression [200, 201]. Previous obser- 
vations of UD patches in co-molding show that matrix tension and tension 
of stitching yarn are the dominating failure mechanisms. Tension failure of 
the UD carbon fibers is hardly observed under processing conditions and 
compression results more likely in instabilities than actual failure. Hence, this 
section only considers damage initialization and progression under tension 
loading. Patches are considered shell-like structures, allowing to apply plane 
stress criteria. 
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6.1.3.1 Damage initialization 


Hashin suggests a smooth piece-wise failure criterion for unidirectional com- 
posites [200, 201]. The initialization of matrix damage in plane stress is 


2 2 
022 012 
—| +1—| =1 , o22>0 (6.6) 
fa Ed 


where 272»? denotes the elastic in-plane transverse Cauchy normal stress and 
012 denotes the elastic in-plane Cauchy shear stress. The tensile strength 
and shear strength of the matrix are given by Y?» and Yj2, respectively. Fiber 
tension damage is initialized at 


2 
O11 
=) =1 , 91ı1>0, (6.7) 
Yu 

where Yj; is the tension strength of the fiber and g4; is the tensile elastic 
Cauchy stress in fiber direction. Carbon fiber damage is neglected due to 
their superior strength. However, the stitching yarn represents also a layer of 
unidirectional fiber reinforcement and its failure is relevant to the character- 
istic failure for UD patches during co-molding. 


6.1.3.2 Damage progression 


Strain-softening damage behavior can lead to mesh-dependent energy re- 
lease in Finite Elements. To alleviate this drawback, the following damage 
progression is formulated in terms of displacement instead of strain utilizing 
the characteristic length l. of an element [202]. The equivalent displacement 
and the initial equivalent stress at onset of matrix tension damage are 


l 
5mo = lc (E22)? -2E?, and omo= DEN ((E22) (022) +2E12012), (6.8) 
MO 


respectively. Here, l denotes a characteristic element length and (e) = 1/2(e+ 
|e|) is the Macaulay bracket, which is used to consider tension stress only. 
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The initial equivalent displacement and the corresponding initial equivalent 
elastic stress at onset of fiber damage are 


lL 
ôro = ley (Ei)? and ofp = = (E211), (6.9) 
FO 


respectively. After damage initialization, the current matrix damage variable 


Dy is defined as 
_ Óut(ÓMt — Omo) 


~ duit(Ome— Ooo) 


Wi 
ôm = le (E52)? +2E2, and öm=2— (6.11) 
O MO 


representing the current equivalent displacement öm: and the final equivalent 


DM (6.10) 


with 


displacement öyr at which the material is completely damaged. The final 
equivalent displacement is computed from the tensile fracture energy of the 
matrix Wy (per area due to the displacement formulation) that is dissipated 
during the fracture. The matrix Young's modulus Ey, shear modulus Gy and 
viscosity 7m are degraded by 


Em=(1-DwWEm , Gu-(ü-Dw)Gw and iw-(-Dw)gw. (6.12) 


Analogously, the evolution of fiber tension damage is described by 


s Örr(ömt — Ógo) 
Ömt(Örr- Öro) 


Wi 
Ort = lo (Eo2)2 and 9-27. (6.14) 


The fiber stiffness Er is degraded by 


D (6.13) 


with 


Ep = (1- Dg)Eg. (6.15) 


This results in the behavior illustrated in Figure 6.2. The equivalent stress 
increases linearly until it reaches 0.9, where e is either F for fiber or M for 
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matrix. Further loading results in damage and linear softening, where the 
final point ó.; is defined by the total fracture energy, i.e. the total area under 
the solid curve. Unloading at partial damage results in the dashed line, where 
the stress reduces with a flatter slope. The damage variable can never be 
reduced, such that a subsequent loading would result in the identical slope. 


6.0 Ont 


Figure 6.2: Relation between equivalent stress and equivalent displacement with linear softening. 


6.1.4 Implementation 


The model described above is implemented in a Simulia Abaqus/Explicit 
subroutine VUMAT. A single ply of patch is represented by two element layers, 
one representing unidirectional carbon fibers and the matrix, and a second 
one representing the stitching yarn in transverse direction. This model allows 
for damage in the matrix due to shear, while the transverse stitching yarn 
remains intact and can be damaged separately due to transverse tension. 
Both layers share all their nodes, as illustrated in Figure 6.3. 


If the damage variable approaches 1, elements may experience significant 
deformation due to the stiffness degradation. Hence, elements are deleted, if 
DM =lor Dr =1. 
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UD layer 
Matrix 


Carbon fiber 


Stitching layer 
No matrix 


Stitching yarn 


Figure 6.3: A single ply of patch is represented by two layers of elements. The fiber orientation 
pr of the UD layer is initially perpendicular to the orientation of stitching yarn. 


6.2 Verification 


The implementation of the patch model is verified for a single element with 
generic parameters to ensure that the model behaves as expected. A single 
layer with a single element 1 mm x 1 mm x 0.1 mm is loaded transverse to 
fiber direction with different fracture energies. The element is loaded beyond 
its strength and then unloaded back to its initial configuration. 


For a transverse stiffness Em = 10GPa and tensile matrix strength Ya = 
500 MPa, results are shown in Figure 6.4. The matrix fracture energy Wy = 
25mJmm ? is twice the nominal elastic energy stored after reaching Y22. 
Therefore the slope after damage initiation should be identical apart from the 
sign. If the fracture energy is reduced to half its value, the damage should oc- 
cur instantaneous after reaching the strength. Increasing the fracture energy 
should result in a flatter slope after reaching the strength limit. The depicted 
behavior in Figure 6.4 agrees with these expectations. Equivalently, correct 
behavior of tension in fiber direction and shear is verified. 
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500 + —— Wm = 12.5m]mm ? 
—— Wm - 25mJmm ? 
400 -| 2 


—— Wn -50mjmm 


Transverse stress 022 in MPa 


0.00 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.10 
Displacement in mm 


Figure 6.4: Single element verification of damage implementation for different matrix fracture 
energies. In all cases, the stress increases with the slope prescribed by the Young’s modulus 
up to the prescribed tensile matrix stress. For Wy = 12.5mJ mm, the energy (area under the 
curve) is just equivalent to the fracture energy at this point and the fracture occurs suddenly. For 
Wu = 25mJmm~, further displacement reduces the force with the same slope, as it previously 
increased. For Wy = 50mJmm_~%, the corresponding slope is flatter, since the energy dissipation 
until complete failure is larger. Unloading results in the expected reduced stiffness for all cases. 
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7.1 SMC 


The primary material under investigation is a glass fiber reinforced SMC 
based on an unsaturated polyester-polyurethane hybrid resin (UPPH) with 
25 mm long glass fibers. The material was developed to improve co-molding 
with unidirectional carbon fiber patches [170]. The employed glass fiber 
multi-end roving has a total Tex of 4800 and consists of 80 strands. Hence, 
cutting one multi-end roving once yields 80 fiber bundles with 200 fibers 
with 14um diameter each, that fall on the SMC carrier foil. The complete 
composition of the material is given in Table 7.1. The material is produced at 
Fraunhofer ICT, Pfinztal, Germany. 


Table 7.1: Composition of UPPH glass fiber Sheet Molding Compound 


Component Trade name Supplier Quantity 
UPPH resin Daron ZW 14141 Aliancys 100.00 parts 
Isocyanate Lupranat M20R BASF 19.50 parts 
Release agent BYK 9085 BYK 2.00 parts 
Peroxide Trignox 117 Akzonobel 1.00 part 
Deaeration aid BYK A-530 BYK 0.50 parts 
Inhibitor pBQ Fraunhofer ICT 0.30 parts 


Glass fiber Multistar 272 Johns Manville 23 vol% (41 wt%) 


It has to be noted that the material under investigation is not subject to a 
rigorous quality control. Hence, the local composition and thickness varies 
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within a single SMC bobbin and especially between different batches. Addi- 
tionally, the properties depend on manufacturing conditions, storage time, 
storage conditions etc. Therefore, the goal of this characterization is to find 
a reasonable approximate description of the material behavior, but not the 
definitive parameterization of the material. 


7.1.1 Thermal properties 


The balance equations for the inner energy of the macroscopic model (4.3) 
and the DBS model (5.6) are closed with the transverse heat conductivity x 
and specific heat capacity cy of SMC. Additionally, the solution requires a 
thermal gap conductance kr to model the heat flow at the contact between 
SMC and mold surfaces. These parameters are identified in this section for 
uncured UPPH-GF SMC material. 


Method. The thermal properties are determined from temperature mea- 
surements in a stack of SMC sheets and were conducted at Fraunhofer ICT, 
Pfinztal, Germany. Ten sheets of 50 mm x 50mm x 1.1 mm are stacked and 
thermocouples are placed centrally between each layer. The first sensor T; is 
placed on the surface of the stack, such that it is positioned between mold 
and stack as soon as the stack comes in contact with the mold. The stack is 
then surrounded with glass wool insulation and placed on a mold surface 
heated to 145°C. A small weight on top of this configuration ensures proper 
contact. The described setup is illustrated in Figure 7.1. 


Results. The measured temperatures T of three samples are shown in Figure 
7.2 as light colored lines for sensors T; to Tg for the first 60s after contact 
initialization. The transient heat transfer is approximated as a 1D process 


according to 
oT | OT 


dr ae 


E (7.1) 


128 


7.1 SMC 


Insulation Ti 


50mm 


Mold heated to 145 *C 


(a) Sketch of the setup (b) Photo of the prepared stack 
Figure 7.1: Setup for the evaluation of the transverse heat conductivity. Ten sheets of SMC are 
stacked with temperature sensors between each layer and placed on a heated plate, while being 
insulated on all other boundaries. 


with ze [0, Ds] and the stack height D, = 11 mm. The specific heat capacity 
cy is approximately 1530] kg! K}, which is estimated from data of the resin 
system [203] and glass by a rule of mixture. The boundary conditions are 


T oT 
mn =-kr(Tm-T) and ET: =0 (7.2) 


zZ 
dz t,z=0 t,z- Ds 


with a constant mold temperature Ty = 145°C. Initially, the temperature is 
homogeneous at T|;=0,z = 24°C. 


An optimal fit to the measured data is obtained by solving the initial value 
problem for an initial guess of x and kr, computing a scalar squared error and 
minimizing the error iteratively. The solution of Equation (7.1) with optimal 
parameters (summarized in Table 7.2) is plotted at the sensor positions in 
Figure 7.2. 


Discussion. The results for senors Ty and Tjo are not considered, as they 
do not show any relevant change during the time frame investigated here. 
The experimental data deviates from the numerical solution of (7.1) due to 
manufacturing related inhomogenities and uneven sheet thickness. The 
fitted parameters are considered an averaged idealization of the heat transfer 
process. The specific heat capacity just affects the heat flux into the SMC, but 
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Figure 7.2: Measured temperatures (light colored lines) and best fit based on one-dimensional 
heat transfer equation (dark colored lines) for sensors Tı to Tg 


Table 7.2: Thermal properties of UPPH-GF SMC in B-staged state. 


Property Symbol Value 


Thermal conductivity K 0.163 Wm! K~} 
Gap conductance kr 403Wm? K! 
Specific heat capacity Cp 1530Jkg ! K! 


not the resulting temperature. Hence, an error in the approximation from the 
rule of mixture does not affect the resulting temperature profiles in this fit 
or the non-isothermal simulations presented later, if the same specific heat 
capacity is used in both cases. The heat flux might be computed inaccurately, 
if the specific heat capacity is inaccurate, but this is of no interest in this work. 


7.1.2 Paste viscosity 


The viscosity of the pure paste (T, y) with additives, but without fibers, is 
of major interest to this work. It is used in the macroscopic reference model 
from Chapter 4 (see Table 4.1) and it is needed in the DBS from Chapter 5. 
For example, the DBS requires the paste viscosity for the matrix contribution 
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in Equation (5.7), to compute drag forces in Equation (5.18) and to compute 
lubrication forces between bundles in Equation (5.25). 


Method. Specimens for viscosity measurement of the SMC paste are prepared 
by filling the paste in a mold instead of processing it on the SMC line after 
mixing. The mold is filled up to approximately 1 mm thickness and sealed 
with a styrene tight foil. The paste is matured for two weeks at room tempera- 
ture similar to the SMC for molding. Round coupons of 25 mm diameter are 
cut from the cured paste and placed in a Anton Paar MCR 501 rheometer at 
Fraunhofer ICT, Pfinztal, Germany in plate-plate configuration. The viscosity 
is measured in oscillatory mode at 20°C, 40°C and 80°C for three specimens 
each. 


Results. The measured viscosity shows typical power-law behavior, as de- 
picted in Figure 7.3. The displayed viscosity refers to the norm of the mea- 
sured complex viscosity. 
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Figure 7.3: Measured viscosity (points) and best fit. 
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As the power-law behavior has a singularity at zero shear rate, a Cross- WLF- 
like model 


T -eiT-T*) 
= MD with MoT) = Diet 13) 


is used to model the paste viscosity. A fixed transition shear rate y = 0.1 is 
employed to limit the viscosity at small shear rates. The power-law coefficient 
n, transition shear rate yo and parameters T*, Dj, aı, a are fitted to the 
experimental results. The fit is obtained by minimizing the squared sum 
of the normalized least-squares error at each measured temperature for all 
specimens. The results of the best fit are illustrated in Figure 7.3 and the 
optimal parameters are summarized in Table 7.3. 


Table 7.3: Viscous properties of UPPH paste in B-staged state. 


Property Symbol Value 


Reference viscosity Di 72kPas 
Transition shear rate Yo 0.1571 
Power-law coefficient n 0.385 
Temperature parameter T* 40.73°C 
Fitting parameter ay 7.94 
Fitting parameter aa 105.96 °C 


Discussion. All computed values of the fit displayed in Figure 7.3 are within 
the measured data range, except for low shear rates at 80 °C. The parameters 
should not be considered strictly Cross-WLF parameters, as the experimental 
data in Figure 7.3 shows no plateau. The model is mainly chosen to smoothly 
transition to a finite viscosity as the shear rate approaches zero. 
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7.1.3 Compaction 


In both models, the macroscopic reference model from Chapter 4 and the DBS 
from Chapter 5, pressure is related to the volumetric strain via an equation 
of state p(Ae). This relation is obtained in this section as tabulated data for 
subsequent simulations. 


Method. SMC stacks of 800 mm x 450 mm size, which represents 10096 mold 
coverage, are consolidated in a Dieffenbacher DYL630/500 hydraulic press 
with parallelism control at Fraunhofer ICT, Pfinztal, Germany. The mold 
temperature is 145°C and the mold closes with a constant speed of 1 mm s^! . 
Automated cutting of SMC sheets ensures an accurate initial mold coverage 
and compression data (force, mold displacements) is recorded during the 
trials. The expected relation between volumetric strain and pressure is illus- 
trated in the schematic in Figure 7.4: Upon first contact, the bulk material 
offers small resistance to compression, as trapped air is compressed and 
released. The resistance increases, as an increasing amount of pores is closed 
and air escapes until the maximum compression force is reached. At this 
constant compression force, the material first expands due to heating and 
consequently shrinks due to cross linking of the thermoset. Finally, the cured 
part expands elastically as the compression force is relaxed and the mold 
opens. The final recorded gap is the part thickness. 


Results. The measurement in Figure 7.5 shows the behavior sketched in 
Figure 7.4 for three samples molded with the identical configuration of three 
SMC sheets each and is plotted with light colors. The thickness of the de- 
molded part is indicated by dots and is in agreement with the measured 
displacement at release. However, the entire mold and press are elastic and 
an empty stroke is used to subtract the pressure-dependent mold displace- 
ment. The so corrected measurements are plotted in dark colors in Figure 7.5. 
Only the region of the rising flank is of interest for the compression molding 
simulation here. Therefore, these points are extracted and plotted over the 
Hencky strain Ae = log (A/ ho) in Figure 7.5. An averaged relation between 
strain and pressure is obtained by averaging the strains at various pressure 


133 


7 Characterization 


Pressure 


Gap 
Figure 7.4: Schematic relation between pressure and mold gap for compression. 1) First contact 


between mold and SMC stack. 2) The maximum compression force is reached. 3) The maximum 
thermal extension is reached. 4) The part is fully cured. 5) The part is demolded. 


levels. This tabulated data (see Table 7.4) is then used to interpolate the 
equation of state p(Ae) in subsequent simulations. In rare cases, where the 
tabulated data range is exceeded, the data is extrapolated linearly. 


Table 7.4: Relation between Hencky strain and pressure for UPPH-GF SMC in B-staged state. 


Hencky strain Ac Pressure p in MPa Hencky strain Ac Pressure p in MPa 


-0.0000 0.0 -0.2029 63.2 
-0.0770 6.3 -0.2073 69.5 
-0.1098 12.6 -0.2116 75.8 
-0.1325 18.9 -0.2167 82.1 
-0.1496 25.3 -0.2219 88.4 
-0.1638 31.6 -0.2270 94.7 
-0.1749 37.9 -0.2317 101.1 
-0.1840 44.2 -0.2349 107.4 
-0.1923 50.5 -0.2378 113.7 


-0.1982 56.8 -0.2407 120.0 
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Figure 7.5: Compressive behavior of the SMC. (a) The recorded data (light color) is corrected with 
the stiffness of the mold and press, which is determined from an empty stroke. The corrected 
data is shown in dark colors and the final part thickness agrees with the recorded mold gap 
measurements. (b) The rising flank is extracted and an interpolation is used as tabulated data to 
describe the relation between Hencky strain and pressure. 


Discussion. The resin system under investigation shows compressible be- 
havior, as previously reported in [155]. The measurement agrees with the 
proposed mechanisms, although the dependence on these specific mech- 
anisms is not explicitly proven. The result may be negatively affected by 
spillage through the mold gap, but the amount of material lost by spillage is 
typically rather small compared to the bulk material. Only the rising flank 
is considered relevant for compression molding simulations, because the 
loading in this application is monotonic. However, the overall behavior is 
irreversible, as seen in Figure 7.5a. Finally, the temperature dependence 
and rate dependence of the paste viscosity may introduce an error to this 
quasi-static analysis, if the ability to close pores is affected by resin flow. 
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7.1.4 Friction at mold surface 


The friction between mold and SMC is an important factor in any SMC com- 
pression molding model and is introduced in Section 4.1.4. It can be deter- 
mined experimentally from pressure differences in a press rheometer, which 
has been introduced in Section 4.2 and Figure 4.4. This section reports pres- 
sure results from molding trials with a press rheometer and discusses several 
friction models in a parametric study with the one-dimensional model from 
Section 4.2. 


Method. The press rheometer is equipped with up to 13 Kistler 6167ASP1.6 
pressure sensors, which can measure pressures up to 200 bar at mold temper- 
atures up to 200°C. The placement of these sensors in the mold is illustrated 
in Figure 7.6 and the outer dimensions of the mold are 800 mm x 450 mm. 


Position y in mm 


0 100 200 300 400 500 600 700 800 
Position x in mm 


Figure 7.6: The pressure sensors P, to Pıo are placed along the center line at 32 mm, 146 mm, 
248 mm, 350 mm, 450 mm, 552 mm, 604 mm, 654 mm, 709 mm and 764 mm. The sensors P1; to 
P13 are positioned in a perpendicular line at 500 mm. 


Trials are performed in collaboration with Sergej Ilinzeer at Fraunhofer ICT, 
Pfinztal, Germany with stack dimensions 600 mm x 450 mm (3 layers, 7596 
mold coverage) and 400 mm x 450 mm (4 layers, 50% mold coverage) on a 


136 


7.1 SMC 


Dieffenbacher COMPRESS PLUS DCP-G 3600/3200 AS hydraulic press with 
parallel cylinder control. The mold is heated to 145°C for all trials and is 
closed with a constant speed h = -1mm s™! until the compression force 
reaches 4400 kN. After reaching the maximum compression force, the mold is 
kept closed for 120s to cure the resin and plates are subsequently demolded. 
The pressures of sensors, displacement of hydraulic cylinders and compres- 
sion force are recorded during the trials. The recorded data of five specimens 
in each configuration is aligned at the switch-over point to force control. This 
point is t = 0 in the evaluation, hence t < 0 means that the press operates in 
displacement control and t > 0 means that the press tries to keep the force to 
a constant value of 4.4 MN. 


The press rheometer is modeled by the one-dimensional model developed in 
Chapter 4, which describes the SMC as an anisotropic, compressible, non- 
Newtonian, non-isothermal material with the parameters determined in the 
previous sections. All parameters of the numerical model except for the fric- 
tion model are determined independently from the press rheometer. Hence, 
the one-dimensional model is used to efficiently vary frictional parameters to 
find a qualitatively and quantitatively fitting model. The objective is to find 
a single parameterization that can describe the effects observed in the 50% 
coverage and 75% coverage trials. 


Results. First, the friction coefficient of a hydrodynamic power-law friction 
model is varied in the range 1.0 MN sm? to 4.0MNsm with a constant 
power-law coefficient m = 0.6 and reference velocity vg = Imms!. These 
are typical values reported for SMC in literature measured after the squish 
phase for stable plug-flows [3, 147, 148, 150]. 


The simulated press data and experimental press data of the 75% mold cover- 
age configuration are shown in Figure 7.7. Regardless of friction parameters, 
the computed final gap is larger than the experimentally observed thickness. 
The initial compression force is similar for all tested friction coefficients, but 
the final time at which the virtual press controller switches to force control is 
earlier for higher hydrodynamic friction coefficients. In the force controlled 
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Figure 7.7: Press data and pressures at 75% mold coverage for different hydrodynamic friction 
coefficients A - experiment (light colors) and 1D model (dark colors). 


phase, the numerical compression force oscillates. The numerical compres- 
sion forces for higher friction are closer to the experimental force and the 
earlier switching point leads to a reduced slope of the gap, which resembles 
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the experiments well, except for the different final thickness. However, com- 
pression force and mold displacement alone do not allow for an unambiguous 
calibration of the friction parameters. The compression force is a holistic 
measure that does not distinguish the contributions from friction and the 
behavior of the SMC itself. Therefore, the pressures at the sensor positions of 
the model are compared to the measured pressures, too. At A = 1.0MNsm^?, 
the spread between the highest pressure and lowest pressure is much smaller 
than in the experiments and the last sensor is triggered too early. An increase 
to A = 2.0MNsm^? increases the spread, but still does not result in a suffi- 


3 and 


ciently large initial slope of the leftmost sensors. For A = 3.0MNsm” 
A 2 4.0MNsm ^, the computed pressures agree well with the experiments in 


terms of spread, slope, absolute pressure level and compression time. 


The results for 5096 mold coverage are depicted in Figure 7.8. The resulting 
gap and part thickness agree and larger friction coefficients result in a dis- 
placement profile closer to the experiment. The compression force increases 
with a similar slope for all tested friction parameters during the initial squish 
phase. Beyond that, the compression force is closer to the experimental 
value for larger values of A, but does not describe the initial hump in the 
force profile. The compression time from first contact to the switch over 
point and complete fill agrees well with the experiment for A = 3.0MNsm ? 
and A = 4.0MNsm^?. The pressure plots corroborate the observations of 
press data. The hydrodynamic friction model is not able to predict the initial 
pressure peak of the experiment with 50% initial mold coverage. However, 
the solution with A = 3.0MNsm ? predicts the slopes correctly in the force 
controlled phase (t > 0). The time at which each pressure sensor is triggered 
agrees with experiments as well, which indicates a correct compressibility 
model. 


The hydrodynamic friction model does not apply any resistance to a stack 
motion relative to the molds at absent relative velocity. The initial hump in 
the experimental data might be caused by such a sticking behavior, though. 
Therefore, also constant friction is implemented in the one-dimensional 
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Figure 7.8: Press data and pressures at 50% mold coverage for different hydrodynamic friction 
coefficients A (Tg = 0 Pa) - experiment (light colors) and 1D model (dark colors). 


model from Section 4.2 (see also Equation (4.49) for the corresponding ana- 
lytical solution) and varied between 30 kPa and 200 kPa. A constant friction 


value of 30 kPa has been previously employed successfully to describe the 


friction of the considered SMC material [6]. The resulting simulated press 
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Figure 7.9: Press data and pressures at 50% mold coverage for different constant friction stresses 


To (A=0N sm) - experiment (light colors) and 1D model (dark colors). 


data is compared in Figure 7.9. The compression force increases rapidly until 
the capability of the sticking friction is exceeded and then increases slower 
with an increase of friction caused only by the increased contact area between 
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SMC and molds. The experimental kink in the compression force curve is 
between the simulated compression forces with 100 kPa and 200 kPa. 


Discussion. At 75% initial mold coverage, the computed final gap is larger 
than experimentally observed gap. As the initial thickness is a fixed value 
and the compressibility is known from Section 7.1.3, this indicates a loss of 
material through the mold gap in the experiment, while the simulated mass 
remains constant throughout the process. The numerical compression force 
oscillates in the force controlled phase, as the solver uses large time steps 
in this phase, which makes it more difficult for the controller to maintain a 
constant force. 


The simulation model with 50% initial mold coverage predicts a monotone 
increase of pressure from the flow front to the other end of the mold due to 
its plug-flow assumption. This is in disagreement with the experiment, where 
the inner layers get squeezed out of the stack and hence lower pressures are 
observed at the first sensor position (green curve) compared to the second 
sensor (orange curve). This will be discussed in more detail in Section 8.3. 


The constant friction allows to modify the compression force at which the 
flow starts, but at high constant friction values, the mold would need much 
more time to fill than experimentally observed. This can be seen by a flat 
displacement slope in the force controlled phase. The pressure sensor results 
show how aconstant friction is able to describe the initial pressure level. How- 
ever, without reduction of this value after the initiation of flow, the pressure 
spread is vastly overestimated. 


This parameter study leads to the conclusion that both, an initial sticking 
behavior with To = 100 kPa to 200 kPa as well as a hydrodynamic friction, 
are required to model the full SMC compression molding process including 
squish flow and a stable plug flow. Hence, the friction model from Equation 
(4.13) is parameterized according to Table 7.5. 
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Table 7.5: Mold friction parameters 


Property Symbol Value 
Stick friction stress TO 150 kPa 
Transition velocity Ut 10mms^! 
Reference velocity vo 1mms 
Power-law coefficient m 0.6 
Hydrodynamic friction coefficient A 3.0MNsm” 


7.2 Patches 


The unidirectional carbon fiber reinforced patches are produced from ZOLTEK 
UD300 stitch bonded non-crimp fabric UD300 with 50K PX35 rovings at 
Fraunhofer ICT, Pfinztal, Germany [170]. Compared to the SMC formulation, 
the composition (see Table 7.6) is designed such that the B-stage is reached 
faster. This enables direct cutting after production of the semi-finished 
patches without a long maturing phase. Like the SMC production, the manu- 
facturing process of semi-finished patch material is not subjected to quality 
control. Local dry spots due to insufficient impregnation, wet spots due to 
insufficient mixing of the components and uneven patch thickness distribu- 
tion occur frequently in the semi-finished material. This is accompanied by 
a dependence on individual batch production conditions and the storage 
time until processing or testing. These conditions lead to large scatter of the 
semi-finished patches. Nonetheless, an attempt is made to approximately 
determine the properties of this material. 


7.2.1 Thermal properties 
Similar to the SMC, transverse heat conductivity x, specific heat capacity 


Cp and thermal gap conductance kr are required to solve the balance of 
inner energy in patches. The parameters are determined analogously to the 
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Table 7.6: Composition of UPPH patches 


Component Trade name Supplier Quantity 
UPPH resin Daron 41 Aliancys 100.00 parts 
Isocyanate Lupranat M20R BASF 25.00 parts 
Impregnation aid BYK 9076 BYK 3.00 parts 
Styrene Mono styrene BASF 2.90 parts 
Release agent BYK 9085 BYK 2.00 parts 
Peroxide Trignox 117 Akzonobel 1.00 part 
Inhibitor pBQ Fraunhofer ICT 0.30 parts 
Accelerator BorchiKat 0248 Borchers 0.17 parts 
Carbon fiber PX35 UD300 ZOLTEK 50 vol% (62 wt%) 


procedure described in the previous Section 7.1 at Fraunhofer ICT, Pfinztal, 
Germany. 


Method. Six patches of 50 mm x 50 mm x 0.45 mm are stacked and thermo- 
couples are placed centrally between each layer. The first sensor Tj is placed 
on the surface of the stack, such that it is positioned between mold and stack 
as soon as the stack comes in contact with the mold. The stack is then sur- 
rounded with glass wool insulation and placed on a mold surface heated to 
145*C. A small weight on top of this configuration is used to lightly press the 
stack against the mold surface. 


Results. The measured temperature profiles are shown in Figure 7.10 with 
light colored lines. The resulting properties of the best fit and the specific 
heat capacity from rule of mixtures are listed in Table 7.7. The solution of the 
one-dimensional heat transfer Equation (7.1) with these optimal parameters 
is plotted in dark colors in Figure 7.10. 


Discussion. The heat transfer from mold to the stack is faster than in SMC, 
but the transverse conduction within the stack is similar to the SMC. A pos- 
sible source of error is the adhesion quality of the patches, which is not as 
good as the adhesion between SMC sheets. Therefore, the transverse heat 
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Figure 7.10: Measured temperatures (light colored lines) and best fit based on one-dimensional 
heat transfer equation (dark colored lines) for sensors T1 to Tg. 


Table 7.7: Thermal properties of UPPH-CF UD patches in B-staged state. 


Property Symbol Value 
Thermal conductivity K 0.128 Wm! K-1 
Gap conductance kr 670W m”? Kl 
Specific heat capacity Cp 1195]kg ! K! 


conductivity could be underestimated due to imperfect interface bonding 
between patch plies. 


7.2.2 Tensile stiffness and strength 


The patch model in Chapter 6 requires the Young's modulus in carbon fiber di- 
rection E; and in stitching yarn direction Er as well as the matrix Young's mod- 
ulus Em. Further, the corresponding strengths Yi v Yi1 and Y22 are needed to 
parameterize the model. Hence, these values are approximated via tensile 
tests at different temperatures in this section. 


Method. Pre-impregnated stitched specimens with dimensions 20 mm x 
160 mm x 0.45 mm are tested under tensile loading with a clamping length of 
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80 mm at KIT IAM-WK, Karlsruhe, Germany. At room temperature, specimens 
are tested on a ZwickRoell Z2.5 universal testing machine with 3 mm min"! 
and 300 mm min“! deformation velocities. For specimens with fiber orien- 
tation in loading direction, testing is started after applying a tension load of 
100 N and the stiffness is evaluated at 0.3% strain. For specimens with fiber 
orientation perpendicular to the loading direction, testing is started after ap- 
plying a tension load of 2.5 N and the stiffness is evaluated at 0.3% strain. Five 
specimens were tested in each configuration, but only four experiments can 
be considered valid for transverse testing, as other specimens experienced 


damage within the clamping area (see Figure 7.11). 


| 302 20% | 


Figure 7.11: Tensile test with fibers oriented perpendicular to the loading direction are con- 
sidered valid if failure occurs outside the clamped region. If progressive damage occurs in the 
clamped region, the test is considered invalid. 


Figure 7.10 shows that the top part of the patch closest to the mold surface 
(orange sensor) reaches a temperature of 70°C to 90°C after 10s. Only the 
sensor in direct contact to the mold exceeds this temperature. As co-molding 
processes are typically shorter than 10s, a temperature of 80?C is considered 
as the upper limit for a homogeneously heated patch for characterization. At 
this temperature, specimens are tested on an INSTRON ElektroPuls E3000 
universal testing machine equipped with a temperature chamber. The same 
deformation rates and initial forces as for room temperature are applied. 
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Results. The results of five specimens per rate tested at room temperature 
and in fiber direction are shown in Figure 7.12a. The strength in fiber direc- 
tion exceeds the capabilities of this machine. The resulting stiffness with 
perpendicular fiber orientations at room temperature is shown in Figure 
7.12b and the strength is shown in Figure 7.12c. 


300 mm min”! | | H I | 
3mmmin” | i H | 
0 5 10 15 20 25 
(a) Tensile stiffness in fiber direction in GPa 
300 mmmin”! | | | | 
3mmmin^ | | | 
0 0.5 1 15 2 2.5 
(b) Tensile stiffness perpendicular to fiber direction in GPa 
300 mm min”! | | | 
3mmmin” | | | | 


0 5 10 15 20 25 
(c) Tensile strength perpendicular to fiber direction in MPa 


Figure 7.12: Tensile stiffness and strength at 20°C. 


The resulting tensile stiffness in fiber direction at 80°C is shown in Figure 
7.13a. Tensile stiffness and strength at 80 °C for fibers oriented perpendicular 
to the load direction are shown in Figure 7.13b and 7.13c, respectively. 


Discussion. The tensile stiffness in fiber direction does not depend on defor- 
mation rate and temperature. It is considered a constant value of Ej, = 14 GPa 
subsequently. The tensile strength in fiber direction is high enough to con- 
sider fracture in the direction of carbon fibers negligible. The tensile stiffness 
perpendicular to fibers shows no statistically significant dependence on the 
deformation rate, but it is temperature dependent. The summed stiffness 
in perpendicular direction is approximately Ep + Ey = 1.3GPa at 20°C and 
0.4 GPa at 80°C. To parameterize the model and separate the yarn stiffness 
from the matrix stiffness, some rough approximations are made: It is as- 
sumed that the transverse stiffness at higher temperatures is dominated by 
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Figure 7.13: Tensile stiffness and strength at 80°C. 


the stitching yarn and hence the stitching yarn stiffness is approximated by 
Er = 0.3 GPa. Similar to the carbon fibers, this fiber stiffness and strength is 
assumed to be independent of temperature, and the strength is conserva- 
tively estimated as 3 MPa. The matrix is assigned a tensile stiffness of 1.0 GPa 
and strength of 16 MPa at 20°C, which are linearly discounted by a factor of 
1/10 for an increase to 80°C. This is certainly just a very rough estimation of 
the tensile properties, but large scatter has been observed between individual 
batches of the patch material depending on the impregnation quality and 
storage conditions. A robust quality control in the production process is 
necessary to determine exact properties. Nonetheless, the obtained values 
give a reasonable estimation to demonstrate the co-molding model. The set 
of parameters in Table 7.8 is applied in later computer simulations for the 
tensile parameterization of the patches. 


7.2.3 Bending stiffness 
Bending in a rheometer enables access to the bending properties, if bend- 


ing stiffness and membrane stiffness are decoupled, as observed in many 
prepregs [199, 204]. This section investigates whether such a decoupling has 
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Table 7.8: Tensile properties used for simulation of carbon fiber UD patch in B-stage 


Property Symbol 20°C 80°C 
Young’s modulus in carbon fiber direction Ej 14GPa 14GPa 
Tensile strength in carbon fiber direction Y 1 "oo" "oo! 
Young's modulus in stitching yarn direction Eg 0.3GPa 0.3 GPa 
Tensile strength in stitching yarn direction Yu 3MPa 3MPa 
Young's modulus of matrix Em 1.0GPa 0.1GPa 
Tensile strength of matrix Yo? 16GPa 1.6GPa 


to be considered for the deformation in the SMC co-molding process, where 
relatively stiff patches are placed initially at room temperature into the mold. 


Method. A photo ofthe bending apparatus in an Anton Paar MCR 501 rheome- 
ter, which is applied for this investigation at Fraunhofer ICT, Pfinztal, Ger- 
many, is shown in Figure 7.14a. The specimen dimensions are 25 mm x 
50 mm x 0.45 mm. They are subjected to nominal rotation rates of 0.1 min", 
1 min”! and 10 min“! at 25°C as well as 0.1 min”! and 10 min"! at 80°C. The 
prescribed deflection angle a ranges from 0? to 60°. The kinematics of the 
apparatus lead to a nominal curvature of the specimen 


Kp = Fo COS Qp + ro9(1 cosa) cot Ap (7.4) 


with the distance from clamps to the pivot point ro = 9.5mm (details may 
be found in Appendix A.1). The Young's modulus for bending may then be 
approximated by beam theory 


Egb = ———— (7.5) 
with the bending moment Mp and the moment of area I; for deflection angles 
Ab #0. 


Results. The results are shown in Figure 7.14b. The initial large values are 
caused by the singularity in Equation (7.5). The subsequent increase of values 
is the build up of tension closing the clearance between the apparatus and 
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specimen, which was larger during the trials at 80°C. The maximum and 
minimum Young’s modulus determined from tensile testing are added as gray 
dashed lines for reference. 


Young's modulus Epp in GPa 


0 10 20 30 40 50 60 
Deflection angle aj in ? 


(a) Photo of the bending setup (b) Results of the bending tests 


Figure 7.14: Rheometer bending test of the patch material. The solid lines in (b) indicate results 
at 25 °C, the dashed lines indicate results at 80°C. 


Discussion. The approximated stiffness at 25°C is even larger than during 
the tensile tests. This can be attributed to either friction in the apparatus, 
or the fact that these specimens come from a different batch. The stiffness 
at 80°C is about 5 MPa, which is rather an underestimation due to the play 
in the apparatus during these tests, which leads to an overestimation of the 
curvature kp in Equation (7.5). In any way, the bending stiffness is in the 
same order of magnitude as the tensile stiffness. The conclusion is that the 
behavior of these patches seems not decoupled and can be described by 
conventional shell elements. 
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7.2.4 Shear stiffness and strength 


In addition to the tensile properties, the shear stiffness Gy, shear strength 
Yı2, shear viscosity nm and fracture energy Wy are required to complete the 
patch model from Chapter 6. These parameters are obtained in this section 
from bias extension tests. 


Method. The shear stiffness and strength are obtained from bias extension 
tests at KIT IAM-WK, Karlsruhe, Germany with a setup similar to that used by 
Schirmaier [205]. Specimens with dimensions 160 mm x 80mm x 0.45mm 
and with 45° fiber orientation are prepared on an automated cutting table 
by Zünd Systemtechnik AG, Altstátten, Switzerland. They are subjected to 
5mm min~, 50 mm min“! and 500 mm min”! deformation rates on a Zwick- 
Roell Z2.5 universal testing machine at room temperature (see Figure 7.15). 
Additional tests are performed at 80°C on a ZwickRoell Zmart.Pro universal 
testing machine equipped with a temperature chamber at 5mm min”! and 
400 mm min“ deformation rates. The free length is 135 mm and adhesive 
tape is employed to improve the fixation in the clamped region. 


The deformation is measured via Digital Image Correlation (DIC) with an 
unstructured pattern of white marker points. The marker points are painted 
on the rovings of the patches to enable an evaluation of the deformation 
between individual rovings (see Figure 7.16a). Contrary to dry fabrics, which 
can be evaluated with a regular grid [205], the points are distributed in an 
unstructured grid, as rovings are deformed from the manufacturing process 
of the patch and not equally spaced. The deformation is recorded with a 
Canon EOS 70D DSLR camera at 30 frames per second with a resolution of 
4 pixels mm^!. The gray-value image is seeded with tracking points based on 
regional brightness maxima of the image (see Figure 7.16b). A MathWorks 
Matlab DIC tool [206] is used to compute the deformation of these tracking 
points with a correlation window size of 21 pixels x 21 pixels. The resulting 
deformation of the tracking points is then imported to Paraview for further 
processing. A Delaunay triangulation is employed to generate a mesh be- 
tween tracked points (see Figure 7.16c). After removing strongly distorted 
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(a) Sketch of the bias extension setup (b) Photo of the bias extension setup 


Figure 7.15: Bias extension test. The region without clamped carbon fibers is subjected to shear 
deformation. 


elements, the displacement can be interpolated with shape functions (see 
Figure 7.16d) and the deformation gradient of each cell is computed. The 
employed deformation measure is the Green-Lagrange strain (see Equation 
(2.80)), which is computed from an element-size weighted average of all 
elements. 


Results. The resulting shear stress at room temperature is shown in Figure 
7.17a for all tested samples. It is plotted over the Green-Lagrange shear strain 
obtained from DIC. 


The scatter between individual specimens in terms of stiffness and especially 
in terms of strength is large. There is a statistically significant rate dependency 
with a steeper initial slope at higher deformation rates (see Figure 7.18a), but 
not distinct rate effect on strength (see Figure 7.18b). The results at 80?C are 
shown in Figure 7.17b for fewer samples. There is also large scatter between 
samples and no distinct rate effect (see Figure 7.19). 
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(a) Marker points (b) Initialization (c) Mesh (d) Evaluation 
Figure 7.16: Steps for DIC analysis. (a) The specimen is prepared with white marker points. (b) 
The tracking points are identified in the initial frame. (c) The resulting deformed points are 
interpolated by Delaunay triangulation (d) The displacement field is interpolated with shape 
functions and deformation measures are evaluated on the evaluation mesh. 
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(a) Results of the shear tests at 20°C (b) Results of the shear tests at 80°C 


Figure 7.17: Shear test results. The deformation of the highlighted dark colored specimen is 
displayed in Figure 7.20. 


The failure at 20°C varies between abrupt fracture in well impregnated 
patches to a progressive damage in patches with lower impregnation quality. 
A representative sample for such a progressive failure is highlighted in dark 
color in Figure 7.17a and the deformation is shown in Figure 7.20. The dis- 
played strain is the maximum principal value of the Green-Lagrange strain 


\/ (Ex - yw" + E obtained from DIC. The specimen fails under shear load 
first. This corresponds to the first peak in Figure 7.17a and can be seen asa 
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Figure 7.18: Shear stiffness and strength at 20°C. 
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(b) Shear strength in GPa 


Figure 7.19: Shear stiffness and strength at 80°C. 


shift between bundles at the localized shear zone in Figure 7.20b. Then the 
stitching yarn fails progressively with the formation of a single shear band 
parallel to the rovings. 


Discussion. The shear stiffness at 5 mm s^! test velocity is 0.110 GPa at 20°C 
and increases apparently with increased deformation rates. For 500 mm s! 
test velocity, the rate Ei: is nominally 0.06 s! and at Ej? = 0.01 the shear 
stress is = 3.5 MPa. This suggests a shear viscosity of approximately 10 MPa s 
with the model in Equation (6.4). The value nm = 10MPas is subsequently 
used as a rough approximation of the matrix viscosity of the patch without 


2 is esti- 


considering shear-thinning effects. The fracture energy of 5 mJ mm^ 
mated from the enclosed area under the stress-displacement curve. Similar to 
the tensile properties, all matrix properties are simply linearly discounted by 


a factor 1/10 for a temperature increase from 20 °C to 80°C. The results may 
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Figure 7.20: Experimental deformation at 20°C for different crosshead travel states. The failure 
causes local loss of some facets, if these cannot be accurately correlated in a picture after cracks 
occur. 


be negatively affected by the introduction of adhesive tape in the clamping 
region, even though the evaluation of strains via DIC negates some of this 
error. The set of parameters in Table 7.9 is used in later simulations for the 
shear parameterization of patches. 


Table 7.9: Shear properties used for simulation of carbon fiber UD patch in B-stage 


Matrix property Symbol 20°C 80°C 

Shear modulus GM 0.110GPa 0.011GPa 
Shear viscosity NM 10MPas! 1MPas™! 
Shear strength Yi? 6MPa 0.6 MPa 
Fracture energy Ww 5mJmm ? 0.5mJmm? 


Simulation. To verify correct qualitative and quantitative behavior of the pa- 
rameterized model, the bias extension test is simulated with the patch model 
described in Chapter 6 and using the material parameters from Table 7.7, Ta- 
ble 7.8 and Table 7.9. The modeling domain is partitioned in clamping regions 
and free region (see Figure 7.21). The free region is further partitioned with 
5mm wide sections representing the width of rovings in the patch. As this 
is a characteristic length of the specimen, the entire domain is meshed with 
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triangular elements of approximately 5mm edge length. The partitioning 
ensures that one edge of each element is always oriented in fiber direction. 


ibed moti 
ene N. horizontally fixed 


horizontally fixed 
Ko 
fixed 


Figure 7.21: Boundary conditions and partitioning of the simulated bias extension test. 


The patch clamping is critical and boundary conditions are chosen carefully. 
The clamping regions are coupled to reference points with kinematic con- 
straints that prohibit any motion in the direction of carbon fiber rovings. How- 
ever, the transverse motion is not constrained, as the experimental clamping 
does not reliably prohibit this motion. Instead, the edges marked at the top 
right and bottom left in Figure 7.21 have constrained horizontal displacement 
to ensure that the upper clamping does not deviate from the loading axis. 
The recorded time-displacement curve from the testing machines are then 
applied to the upper reference point, while the lower reference point remains 
in a fixed position. 


The resulting deformation behavior is illustrated in Figure 7.22. First, the 
shear band forms and leads to a localized shift between bundles in Figure 
7.22b. The position of the shear band is determined by numerical or experi- 
mental imperfections. Its formation is observed either at the upper position 
or at the lower position. Even after complete failure of the matrix material, 
the stitching yarn still holds both halves of the specimen together (see Figure 
7.22c). Finally, the stitching breaks completely and the specimen is sepa- 
rated (see Figure 7.22d). This behavior qualitatively matches the observed 
experimental behavior shown in Figure 7.20. 
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(d) 7.1mm 


(c) 4.7mm 


(b) 42mm 
Simulated deformation of a patch under bias extension loading at 20°C and 


(a) 2.9mm 


Figure 7.22 


The colormap shows the maximum principal value of the Green-Lagrange 


500 mm min-!. 


strain with pink color representing values that exceed the color range. 


The corresponding force-displacement curves are shown in Figure 7.23, 


where the extracted positions of Figure 7.22 are marked by dots. The simula- 


tion agrees with those experimental results showing a progressive damage 


and the images at the extracted positions show that the same mechanisms 


force and damage energy are close to the 


aximum 


apply. The initial slope, m 


ts. Evaluations at 80?C and different rates show similar results 


experimen 


and it is concluded that the model described in Chapter 6 and the calibra- 


tion of parameters from this chapter allow a reasonable prediction of patch 


deformation and damage initialization. 
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Figure 7.23: Force-displacement curves of the bias extension test at 20°C and 500 mm min“! 
Simulation results with the parameterization of this chapter are plotted on top of the experi- 


mental results (light colors) in dark colors and the points indicate the position of the snapshots 
shown in Figure 7.22. 
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In contrast to the previous verification in Chapter 5 and Chapter 6, the vali- 
dation compares the model predictions to experimental results. Therefore, 
parts are molded with a wide range of geometries to explore merits and weak- 
nesses of the proposed mesoscale modeling approach. An overview over the 


validation cases is given in Table 8.1. 


Table 8.1: Overview over the chosen validation cases. 


Simple Honeycomb Press Complex IRTG demon- 
Dome ribs rheometer dome strator 
Section 8.1 8.2 8.3 8.4 8.5 
SMC material UPPH-GF VE-GF UPPH-GF UPPH-GF UPPH-GF 
Patch material - - - UPPH-CF - 
Length 120mm 310mm 800mm 120mm 800mm 
Width 95mm 310mm 50mm 94mm 250mm 
Bundle length 25mm 50mm 25mm 25mm 25mm 
Bundle count 715,000 ~17,000 ~40,000 ~18,000 ~150,000 
Paste viscosity Newtonian Newtonian Cross-WLF Cross-WLF Cross-WLF 
Compressible x x v v v 
Non-isothermal x x v v v 
Evaluation Flow front Flow front Press force Flow front Flow front 
Orientation FVC Pressures FMS FMS 
FMS Press force Orientation Co-molding Orientation 
Curvature Press force 
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8.1 Double curved dome with isothermal 
Newtonian matrix! 


The first validation example should show the general ability of DBS presented 
in Chapter 5 to predict the fiber architecture in a 3D geometry under simpli- 
fied Newtonian isothermal conditions. Therefore, a small geometry is molded 
from SMC with different initial stack configurations and the resulting parts 
are scanned by computer tomography. The real fiber architecture is then 
compared to DBS and to the results of a simplified version of the macroscopic 
model presented 4. 


8.1.1 Experimental setup 


The material under investigation is the glass fiber reinforced UPPH SMC 
that has been described in Section 7.1. The specimen under investigation 
is a curved dome with outer dimensions 120 mm x 94 mm and a nominal 
thickness of 2 mm. The specimen is molded from two different initial stack 
configurations. The split configuration is molded from two stacks with di- 
mensions 80 mm x 30mm x 5.3 mm, which are placed at both ends of the 
mold, as illustrated in Figure 8.1a. This configuration leads to a knit line in 
the center, where both stacks meet. The bundle architecture in vicinity of 
this knit line is of particular interest, as it would represent a weak spot in a 
structural application. The asymmetric configuration is molded from a single 
stack with dimensions 80 mm x 60 mm x 5.3 mm, which is placed on one side 
ofthe mold and enables a longer flow path, as shown in Figure 8.1b. The mold 
is heated to 145 °C and closed with a Lauffer hydraulic press in collaboration 
with Lucas Bretz at KIT wbk, Karlsruhe, Germany. The maximum press force 
was limited to 50 kN. 


] Parts of this section are based on [195]. 
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(a) Split stack configuration (b) Asymmetric configuration 
Figure 8.1: The molded part has outer dimensions 120 mm x 94mm. For the split stack configu- 
ration, two SMC stacks are placed at both ends of the mold. For the asymmetric configuration, a 
single stack is placed on one side of the mold. 


The molded samples were analyzed by volumetric imaging using an Yxlon 
X-ray UCT system with a Perkin Elmer flat panel Y.XRD1620 detector and a 
reflection tube by Comet. A volumetric gray scale image is reconstructed 
based on the Feldkamp, Davis and Kress (FDK) algorithm [207] with a result- 
ing voxel size of 69 um. UCT scans and volumetric image reconstruction were 
performed by Ludwig Schóttl at KIT IAM-WK, Karlsruhe, Germany. 


8.1.2 Numerical setup 


The experimental process is modeled with DBS and a macroscopic model 
according to Chapters 4 and 5. The molding domain is represented by Eu- 
lerian elements and only those Eulerian elements occupied by initial stack 
positions are initially filled with material. 


Both mold halves are represented by rigid shell elements. They interact with 
the SMC paste through normal contacts and hydrodynamic friction (see 
Equation (4.11)). While the lower mold is constrained at a fixed position, the 
upper mold is closed with the profiles given in Figure 8.2. These profiles are 
an idealization of experimentally recorded press profiles. The variation in the 
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experimental profiles can be attributed partly to a non-uniform thickness of 
SMC sheets and to the reaction time of the press control unit. The simulation 
stops after a complete fill with the final part height and does not include the 
subsequent holding and curing process. 


20 X 
Split configuration (experimental) 
——— Split configuration (idealized) 
E 15- Asymmetric configuration (experimental) 
3 Asymmetric configuration (idealized) 
= 
a 104 
S 
dp 
= 
© 
= 54 
0 i i i i i 
0 1 2 3 4 5 6 
Time tin s 
Figure 8.2: Distance between upper and lower mold during the compression flow phase of SMC. 


Six parts of the split configuration were produced and are shown with light green lines. Four 
parts of the asymmetric configuration were produced and are shown with light orange lines. 
Additionally, the idealized mold profiles for simulations are shown in solid green and solid orange 
for split stack configuration and asymmetric configuration, respectively. [195] 


The macroscopic reference solution was obtained with a simplified (New- 
tonian, no anisotropic interaction terms) version of the model described in 
Chapter 4 that resembles the state of technology. The initial fiber orientation 
is described by a planar isotropic fiber orientation tensor. 


For the DBS, bundles are generated with the procedure described in Section 
5.1.2, assuming a planar isotropic orientation distribution in the stacks. The 
total number of bundles per stack is 7600, which is computed from the bundle 
properties and a nominal bundle volume in the part of 5844 mm? at the given 
nominal fiber volume fraction. Each bundle is discretized with ten linear 
truss elements. 
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The explicit time integration requires a small time increment due to the high 
resin viscosity. The mass of the entire model was therefore scaled by a factor 
Km to improve the time increment, while ensuring that kinetic energy remains 
negligible small compared to the external work. A summary of the simulation 
parameters is listed in Table A.6. 


8.1.3 Results and discussion 
8.1.3.1 Flow front progression 


Figure 8.3 provides an overview on the compression molding process simula- 
tion of the split stack configuration. The initial mold gap at t = 0s is 20mm 
and the lowest part of the upper mold is just not touching the SMC stacks. 
Closing the mold initially deforms the stacks into the curved geometry, but 
does not start material flow. During forming, the two-way coupled approach 
pulls the stack sideways in the dome-shaped mold. This can be observed by 
the lateral deformation of the stack tips depicted at t 2 2s in Figure 8.3. The 
mold gap is reduced to the initial stack height of 5.3 mm after approximately 
two seconds. From there on, flow dominates the mold filling process and 
fiber bundles are carried with the SMC until the final part thickness of 2 mm 
is reached. 


8.1.3.2 Orientation and separation effects 


Figure 8.4 shows slices through the midplane of the upper and lower planar 
regions of the scanned part in split stack configuration. The white strands 
represent fiber bundles, which remain in their bundled structure even for 
the applied high degree of deformation. The knit line features a severe fiber- 
matrix separation and only a small amount of fiber bundles bridges the gap in 
this zone. The inner slice in Figure 8.4 even shows some pores. Regions close 
to the mold boundaries and the knit line show a bundle alignment parallel to 
the boundary. Bundles perpendicular to the boundary are likely pulled out of 
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Figure 8.3: Snapshots of the molding process for the split stack configuration. The compression 
molding process starts with a deformation of the two initial stacks. Subsequently, the SMC is 
forced to flow until the part reaches its final thickness of 2 mm. The DBS approach lets bundles 
deform and flow with the matrix material while enforcing two-way coupling. Therefore, the flow 
is naturally anisotropic and depends on the current bundle configuration. [195] 


this region by forces acting over the entire length of the bundle and parallel 
bundles remain close to the boundaries (see also Figure 2.9). Regions farther 
away from boundaries show a regular random in-plane orientation. 


The DBS result is sliced in the same planes and the result is depicted in 
Figure 8.5. The simulation results show a slightly larger area of fiber-matrix 
separation and bundles do not bridge the resin-rich knit line. This might be 
caused either by the experimental setup, because the part was compressed 
further than the nominal thickness, or by an underestimation of the drag 
forces on bundles in the model. Similar to the uCT-scan, boundary regions 
show a predominant orientation parallel to the boundaries. 
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Figure 8.4: Slices through the upper and lower planar regions of the uCT Scan. Fiber bundles 
stay intact during molding and fiber-matrix separation can be observed at the knit line. The knit 
line region includes pores marked with red circles. [195] 


For a quantitative comparison of the DBS to the macroscopic simulation and 
the uCT scans, bundle orientations are evaluated on a uniform 12 x 16 x4 grid 
of sub-volumes. The discrete second-order fiber orientation tensor for each 
of the sub-volumes is computed with the procedure described in Section 
5.1.6.3. The slices of the uCT scan shown in Figure 8.4 are analyzed in 2D 
using the Orientation] plugin of the image processing software Fiji [208]. This 
procedure assigns a major direction to each 10 x 10 pixel area. The discrete 
fiber orientation tensor is evaluated on an equivalent 12 x 16 grid to represent 
the orientation state as tensor components. 


A comparison of the DBS approach, uCT scan and the conventional macro 
fiber orientation model is depicted in Figure 8.6 for the split stack configu- 
ration. The A,,-component of the uCT-analysis features three significantly 
higher oriented vertical stripes at both ends of the mold and the knit line. 
Conversely, the Ay,-component of the uCT-analysis indicates a dominant 
orientation in horizontal direction at the top and bottom mold boundaries 
with lower values at the vertical mold boundaries to the left and right of the 
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wi 


Figure 8.5: Slices through the planar regions of the DBS result. Each gray cylinder represents 
a bundle segment consisting of 200 individual fibers. The knit line at the center is matrix rich 
and no bundles gap this region. Bundles close to the boundaries show a reduced fiber volume 
fraction and more bundles are oriented parallel to the boundary. [195] 


figure. The corresponding DBS is able to reproduce these three stripes of 
higher vertical orientation at the correct positions. Characteristic gradients 
and the level of orientation is predicted well. The macroscopic model using 
fiber orientation tensors and Jeffery’s equation does not account for the con- 
straints at mold walls and shows a homogeneous orientation distribution. 
In homogeneous regions, such as the inner slice with some distance to the 
knit line, Jeffery's equation leads to a reasonable prediction of the orientation 
state. 


The DBS limits any bundle orientation normal to the molds, because bun- 
dle segments cannot be physically arranged in normal direction in the con- 
strained mold gap. Thus, the Az; -component is small in the planar regions of 
the part. An investigation of a magnified uCT scan with higher resolution con- 
firms that fiber bundles at the knit line are primarily oriented in-plane. The 
computation based on fiber orientation tensors shows a dominant normal 
component of fiber orientation at the knit line. 
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Figure 8.6: Comparison of DBS results with uCT Analysis and fiber orientation tensor based 
computation utilizing Jeffery's equation for the split stack configuration. The first row shows 
orientation tensor component Axx which indicates vertical fiber orientation in this representa- 
tion. The second row shows the Ayy-component representing horizontal fiber orientation. The 
third row shows the A7;-component representing fiber orientation normal to the observation 
plane. The orientation analysis of the uCT image slices is limited to two dimensions. Thus, the 
central image in the third row shows a high resolution uCT scan of the region indicated in the 
illustration above. The magnified view reveals a dominant in-plane orientation of bundles. [195] 


Figure 8.7 is analogous to Figure 8.6, but describes the evaluation of the 
asymmetric stack configuration with a maximum flow path of 60 mm in y- 
direction. This configuration confirms observations of the previous case with 
significantly higher orientations parallel to mold walls that are not described 
by tensor based theory. Despite a longer flow path, the magnitude of re- 
orientation is similar to the split stack configuration due to a similar stretch 
in y-direction (50% initial mold coverage each). 
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Figure 8.7: Comparison of DBS results with uCT Analysis and fiber orientation tensor based 
computation utilizing Jeffery's equation for the asymmetric stack configuration. Refer to Figure 
8.6 for a detailed explanation of the layout. [195] 


8.1.3.3 Bundle curvature 


The curvature of bundles is evaluated with the procedure described in Section 
5.1.6.1 at each node k. A color-coded scatter plot of the curvature for the split 
stack configuration is plotted in Figure 8.8. It shows that the largest curvatures 
occur at corners and close to the knit line. The curvature at the knit line 
originates probably from a flow in x direction compressing bundles to a zigzag 
shape. The curvature in the uCT scan is obtained only for the central region 
in order to have sufficient resolution for tracking bundle curvature [209]. 
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Figure 8.8: Simulation results of bundle curvature. The highest values occur at the corners of 
the mold and at the knit line. The part's three-dimensional shape is visible in this plot due to 
the bending of bundles at curvatures of the geometry. High-resolution uCT curvature data for 
the central area marked with a black rectangle is obtained by Ludwig Schóttl at KIT IAM-WK, 
Karlsruhe, Germany. [195] 


The projection of curvature values in the uCT scan region on the y direction 
is plotted in Figure 8.9. The maximal values of the uCT scan agree well with 
the maximal curvatures computed from the DBS. The mean curvature of DBS 
slightly underestimates the curvature in the outer regions, but agrees well 
around the knit line. It Should be mentioned that simulated curvature might 
depend on the segment length of bundles. 


8.1.3.4 Number of contacts 
An a priori estimate for the number of contacts per bundle segment is given 
by equation (2.51). This estimate predicts about 6.2 contacts per bundle 


segment. The total number of contacts (g « 0 in Equation (5.23)) in the 
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Figure 8.9: Curvatures projected to the y axis in the scanned central region (see rectangular 
section in Figure 8.8) encompassing the knit line. 


DBS is evaluated for each frame of the simulation results and is plotted in 
Figure 8.10. This averages to approximately 1.8 x 10? contact pairs for the 
split stack configuration and 2.2 x 10? contact pairs for the asymmetric flow, 
which has a slightly increased fiber volume fraction compared to the nominal 
value. Considering the total amount of 77,438 and 87,950 bundle segments, 
this evaluates to 4.6 and 5.0 contacts per bundle segment, respectively. This 
evaluation fits well with the estimate given in equation (2.51). 


8.1.4 Concluding remarks 


The mold filling process of a double curved dome structure is described 
with DBS using an isothermal Newtonian paste and an isotropic Newtonian 
reference model. The predicted fiber architecture agrees well with measured 
„CT data w.r.t. fiber miss-alignment at mold walls, knit line formation and 
curvature. The flow of material is assumed to be isothermal in this validation 
case. This assumption is quite common for the bulk material of SMC, as 
the time scale of thermal diffusivity in SMC is large compared to the time 
it takes the material to flow (less than 5s). However, even small changes 
in temperature can cause a relevant change in viscosity and is therefore 
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Figure 8.10: Number of bundle-bundle contact pairs during the molding process. The number 
of contact pairs decreases during the forming phase of the stack and increases during flow, when 
the entire stack is compressed. [195] 


introduced to later validation examples. Additionally, lubrication between 
bundles is neglected in this validation example for computational efficiency, 
since it has little effect on fiber architecture (see Section 5.2.3). The DBS is 
able to predict fiber-matrix separation effects at the knit line and thus enables 
a better description of structural weak spots in such areas. At regions close to 
the mold walls and the knit line, DBS accounts for spatial constraints of the 
fiber orientation due to mold boundaries and leads to more accurate fiber 
orientation results than orientation tensor based macroscopic simulation 
models. Nonetheless, Jeffery's equation leads to reasonable results in planar, 
homogeneous regions at up to ten times faster computational times. 


8.2 Honeycomb ribs with isothermal 
Newtonian matrix 


Fiber-matrix separation is investigated in a plate mold with honeycomb 
shaped ribs. The experiments were conducted as part of the Master's thesis 
by Florian Rothenhäusler at Volkswagen AG, Wolfsburg, Germany using a VE 
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SMC with 50 mm fiber length [63].? This section reports the application of the 
DBS model to this geometry and material system. The resulting fiber-matrix 
separation and compression forces are compared to the experimental results. 


8.2.4 Experimental setup 


The material in this example differs from the default material described in 
Section 7.1. Here, the commercial SMC HUP EJ 43529 by Polynt Composites 
Germany GmbH is used for compression molding. It has a fiber volume 
fraction of 36% and fibers are 50 mm long with a fiber diameter of 20 um. The 
mold is a plate with four insert positions, of which one is filled with a rib 
structure consisting of thick ribs (2.5 mm width at top, 2? draft angle), one 
is filled with a thin rib structure (2 mm width at top, 1? draft angle) and the 
remaining two insert positions are filled with dummy plugs. The SMC is 
molded in the two different configurations shown in Figure 8.11 on an RSP 
160.42SZ hydraulic press by Rócher GmbH&Co KG Maschinenbau, Germany. 
In configuration A, the stack is placed directly under the honeycomb ribs and 
the ribs are filled at the beginning of the compression. In configuration B, the 
stack is placed opposite to the honeycomb ribs and the material flows along 
the mold before filling the ribs. Both configurations use a mold closing speed 
of 10 mm s™}, mold temperature of 140°C and maximum compression force 
of 1500 kN. Press forces, velocity and displacement are recorded during trials. 


The fiber volume fraction in the ribs is analyzed by thermal gravimetric 
analysis (TGA). To limit the experimental effort, only the ribs in the position 
marked with labels 1 to 6 in Figure 8.11 are analyzed. This position is chosen, 
because these ribs show a challenging fiber-matrix separation pattern, but are 
also reliably filled during the process. Each of the six walls is subdivided in a 
lower, central and upper partition and TGA is performed following procedure 


2 Parts of the Master's thesis by Florian Rothenhäusler are included here for comparative 
purposes. They will be published soon under the title Experimental and Numerical Analysis 
of SMC Compression Molding in Confined Regions — A Comparison of Simulation Approaches. 
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(a) Configuration A (b) Configuration B 
Figure 8.11: Mold configurations of the honeycomb experiments. 


B of DIN EN ISO 1172:1998. The specimen is pyrolyzed for 30 min at 650 °C 
and subsequently the mineral filler is dissolved by hydrochloric acid [63]. 
This procedure is repeated three times for each configuration (324 samples 
in total) and the fiber volume content is calculated from the fiber mass after 


the extraction procedure. 


In addition to the fully molded samples, short shot moldings in configuration 
A were produced by Simon Wehler at Volkswagen AG, Wolfsburg, Germany. 
The short shots were obtained by inserting 1.6 mm and 4.8 mm thick rigid 
spacers in the mold gap prohibiting complete closure ofthe mold. The surface 
of the resulting parts were scanned with a 3D optical measurement system 
ATOS by GOM GmbH, Braunschweig, Germany. 


8.2.2 Numerical setup 


The process is modeled with the Direct Bundle Simulation technique ex- 
plained in Chapter 5. The upper and lower mold surfaces are extracted from 
the CAD model of the part, placed at an initial distance of 12 mm and are 
meshed with rigid body elements of 2 mm average edge length. The lower 
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mold remains fixed, while the upper mold is either controlled by a virtual 
press controller (see 4.1.5) following the press profile (configuration B), or 
subjected to a constant force (configuration A). Configuration A is simulated 
with a constant compression force, because the press comes effectively to a 
stop after reaching the stack and proceeds in a force controlled mode in the 
trials with this configuration. This is likely caused by the ribs clinging to the 
SMC stack, which delays the start of the flow process. The Eulerian domain is 
expanded upward and downward beyond the mold surfaces and meshed with 
1mm x 1mm x 1mm EC3D8R elements in the region of the honeycombs and 
1mm x2mmx2mm else. The in-plane boundaries are not implemented via 
rigid body contact, but by constraining the velocity in normal direction to 
zero. The regions of initial stacks are filled with material in the beginning and 
17000 bundles of 50 mm are generated in the stack to represent the desired 
fiber volume fraction. The bundle cross section area is set to 1.768 x 1077 m? 
assuming 560 fibers of 20 um thickness each. 


The simulation is isothermal and quasi-incompressible due to a lack of de- 
tailed material characterization of the material employed in this example. 
The paste viscosity is set to a constant Newtonian value of 25 kPas, which was 
obtained as the zero-shear viscosity for the paste employed in this SMC and 
applicable for shear rates up to approximately 10s! [63]. Bundles interact 
with each other, but tangential lubrication forces are disabled as the resulting 
time step would become prohibitive. The friction between mold and SMC 
is modeled with a conventional hydrodynamic friction model that has been 
parameterized with a press rheometer at Fraunhofer ICT, Pfinztal, Germany 
following the procedure described by [150]. A summary of all simulation 
parameters can be found in Table A.7. 


Both configurations are simulated three times each, with different initial 
random fiber bundle distributions. This allows to investigate the effect of the 
initial configuration on process results and allows to quantify the scatter of 
computed fiber volume fraction. 
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8.2.3 Results and discussion 
8.2.3.1 Flow front progression 


The snapshots of the flow front progression in Figure 8.12 show the rib filling 
process for both configurations. 


t=0s 


(a) Configuration A (b) Configuration B 
Figure 8.12: Snapshots of the honeycomb mold filling process. 


In configuration A, ribs are filled shortly after the start of the compression 
process. There is a slight tendency for ribs down the flow path to fill later, but 
the overall filling process occurs quite homogeneous. After 0.8 s simulation 
time, the ribs are completely filled and the flow front in the plate region flows 
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until it reaches the end of the mold. Long fibers, which are entangled in the 
ribs, influence the flow and lead to a significant alignment of fiber bundles in 
flow direction. This can also be seen in Figure 8.13, where the matrix of an 
entire plate was burned off and is compared to the simulated fiber bundle 
architecture. Both, the simulation and pyrolyzed part, feature a distinct fiber 
bundle alignment left to the ribs. 


(a) Pyrolysis [63] (b) Simulation 
Figure 8.13: Comparison of the fiber architecture in configuration A. 


In configuration B, the stack elongates first similar to the linear compression 
flow in a press rheometer. As soon as the flow front reaches the ribbed section, 
a portion of the flow diverges in the ribs and the ribs closer to the initial stack 
are completely filled before the ribs further away even begin to fill. The 
diverged flow causes the formation of a knit line behind the ribs, which can 
be also seen in the bundle architecture in Figure 8.14. 


To validate the flow front progression, the filled domain is compared to 3D 
optically measured surfaces of short shots of configuration A, which are 
shown in Figure 8.15. Figure 8.15a is obtained from a specimen with 4.8 mm 
thick spacers in the mold gap stopping the compression flow soon after 
the stack is touched. Overlaying the simulation time step closest to this 
thickness (see Figure 8.16a) shows that the flow front progression is captured 
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(a) Pyrolysis [63] (b) Simulation 
Figure 8.14: Comparison of the fiber architecture of the thick honeycomb in configuration B. 


accurately. For a later compression state obtained with a 1.6 mm thick blank, 
the simulation and scan also agree (see Figure 8.16b). The simulation predicts 
equal flow front progression in both, thin and thick, ribs. However, the thinner 
ribs appear to be filled slightly faster in the experiment. This is likely caused 
by higher temperatures in the thin region that result in lower viscosity and a 
faster filling progress, which is not accounted for in this isothermal simulation. 
Further, the filling process may be affected by missing parallelism control of 
the press, or the compressibility of the stack. 


(a) Short shot 4.8 mm (b) Short short 1.6 mm 
Figure 8.15: Surface scans of short shots obtained by Simon Wehler from Volkswagen AG, Wolfs- 
burg, Germany. 
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(a) Short shot overlay 4.8 mm (b) Short short overlay 1.6 mm 
Figure 8.16: Comparison of scanned short shots (gray surface) with paste fill region (green 
volume) of thick ribs in configuration A. 


8.2.3.2 Fiber volume fraction 


The focus of this evaluation are the thick honeycomb ribs labeled with num- 
bers in Figure 8.11. These are regularly filled and the observations are qualita- 
tively equivalent to the mirrored ribs in the thin honeycomb. An evaluation 
mesh representing the experimental partitions is used to compute fiber vol- 
ume fractions of the DBS (see Section 5.1.6.3). 


A comparison between thermal gravimetric analysis and the simulation re- 
sults for configuration A is shown in Figure 8.17. The numbers at the bottom 
indicate the rib label according to Figure 8.11 and the + values indicate the 
standard deviation between samples at each position. The experimental 
results show a homogeneous fiber volume fraction distribution in the ribs 
with small deviation between the three analyzed samples. The ribs, which are 
placed directly above the initial stack position (1,2,3) exhibit slightly higher 
fiber volume fractions than the other ones. Contrary, the corresponding 
simulation predicts severe fiber-matrix separation with local fiber volume 
fractions well above the experimentally observed values. An investigation of 
the simulation progress reveals that the ribs fill homogeneously first, but fiber 
bundles are pulled out subsequently. Fiber bundles entangle at the bottom of 
ribs with multiple adjacent ribs (e.g. 1,3,6) and lead to a concentration at the 
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(b) Fiber volume fraction obtained by Direct Bundle Simulation 
Figure 8.17: Fiber volume fraction distribution in thick ribs of configuration A. 


bottom of these ribs. The discrepancy between simulation and experiment is 
likely caused by missing shear thinning properties of the matrix material in 
the model, which lead to smaller forces on the bundle ends extending to the 
base plate. Another reason could be the gap correction described in Section 
5.1.5.3, which is just a coarse approximation of the actual bundle shape. 


The experimental result of configuration B predicts a gradient in volume 
fraction from the top left to the bottom right in Figure 8.18a. The top sections 
of ribs 1 to 3 exhibit large experimental uncertainties - up to 13.3% standard 
deviation at a mean value of 16.6%. The fiber architecture in these regions 
depends strongly on individual realizations of randomness in the mold filling 
process. The numerical model predicts overall slightly higher fiber volume 
fractions and severe separation in Rib 2, that occurs reliably in all three 
repetitions. 


179 


8 Validation cases 


100% 
75% 
50% 
25% 
0% 


100% 
75% 
50% 
25% 
0% 


(b) Fiber volume fraction obtained by Direct Bundle Simulation 
Figure 8.18: Fiber volume fraction distribution in thick ribs of configuration B. 


8.2.3.3 Compression forces 


The simulated compression force, velocity and mold gap are compared to 
press data in order to evaluate the model's ability to predict not just the 
fiber architecture, but also resulting press parameters. The press used in the 
experiments is not designed for precise measurements, but for industrial 
manufacturing processes, and does not have parallelism support. For exam- 
ple, the closing velocity differs from the prescribed value of 10 mm s^! during 
the displacement control and the compression force is not reliably held at 
the nominal value of 1500 kN. This can be seen from the gray experimental 
recordings shown in Figure 8.19 for configuration A and in Figure 8.20 for 
configuration B. The dots in both figures indicate the snapshots positions 
for an easier comparison to Figure 8.12. In configuration A, the press effec- 
tively stops at contact with the stack and tries to proceed force controlled. 
Hence, the DBS model also utilizes a constant compression force in this case. 
This results in a rapid initial compression similar to the experimental veloc- 
ity. However, the subsequent filling phase is faster than the experimental 


180 


8.2 Honeycomb ribs with isothermal Newtonian matrix 


records indicate, which is partially caused by the lower than nominal com- 
pression force in the experiment and partially caused by skipping the initial 
pre-compression phase. In configuration B, the virtual press controller fol- 
lows the averaged press velocity profile and switches to force control, after 
the maximum force of 1500 kN is reached. The resulting displacement dur- 
ing the force control phase is comparable to the experimentally observed 


displacements. 
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Figure 8.19: Compression press data in comparison to DBS for configuration A. Dots mark the 
snapshot positions shown in Figure 8.12. 
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Figure 8.20: Compression press data in comparison to DBS for configuration B. Dots mark the 
snapshot positions shown in Figure 8.12. 
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8.2.4 Concluding remarks 


The molding of a plate with honeycomb shaped ribs has been simulated 
with DBS employing isothermal, Newtonian paste viscosity. The flow front 
progression is reasonable and verified by short shorts. The thin honeycomb 
is filled slightly faster, which can be likely attributed to the neglected tempera- 
ture change. The predicted fiber volume fraction in configuration A disagrees 
with the experimental result, as the effect of fiber bundle pull out is too large 
compared to the experiments. Conversely, the fiber volume prediction in 
configuration B is able to correctly identify critical regions at the top of Rib 1 
and Rib 2. A comparison to other simulation tools (Moldex, 3DTimon DFS) 
showed that macroscopic models are not able to predict the fiber-matrix 
separation in this example and that a neglection of fiber-fiber interactions 
can lead to unreasonably high fiber volume fractions, because multiple fibers 
can occupy the same space [63]. The consideration of non-isothermal con- 
ditions and a non-Newtonian paste seems required to accurately predict 
compression forces as well as bundle architectures in thin structures. 


8.3 Press rheometer 


A press rheometer is a plate mold for rheology experiments at component 
scale and has been introduced in Section 4.2, Figure 4.4 and the characteriza- 
tion of friction in Section 7.1.4. Typically, such a mold is used to characterize 
the macroscopic elongational viscosity and hydrodynamic friction properties 
of the stable plug-flow with an array of pressure sensors that is distributed 
within the mold. The macroscopic elongational viscosity can be measured 
only in such a large mold because only this way it can be ensured that the 
fiber length is sufficiently short compared to the dimensions of the measure- 
ment instrument. Conversely, the proposed DBS model only requires the 
paste rheology and thus the press rheometer is used in terms of validation to 
investigate whether the observed pressures agree with the bottom-up DBS 


182 


8.3 Press rheometer 


model. This is done with precise parallelism controlled presses and the full 
non-isothermal, non-Newtonian, compressible material model. 


8.3.1 Experimental setup 


The experimental setup is already described in Section 7.1.4 and the experi- 
ments are extended by additional trials with the stack configuration 200 mm 
x 450 mm (8 layers, 25% mold coverage) on a Dieffenbacher DYL630/500 
hydraulic press with parallel cylinder control in collaboration with Sergej 
Ilinzeer at Fraunhofer ICT, Pfinztal, Germany. The part and the initial stack 
configurations are sketched in Figure 8.21. The pressures of sensors, displace- 
ment of hydraulic cylinders and compression force are recorded during the 
trials. The recorded data of five specimens in each configuration is aligned at 
the switch-over point to force control. This point is t = 0 in the evaluation, 
hence ź < 0 means that the press operates in displacement control and t > 0 
means that the press tries to keep the force to a constant value of 4.4 MN. 


800 mm 


uiu Qcr 


Figure 8.21: Mold configurations of the press rheometer experiments. 
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8.3.2 Numerical setup 


The SMC flow is dominated by a horizontal one-dimensional elongation and 
the variance in transverse direction is expected to be rather small. Hence, for 
computational efficiency, the DBS domain width is reduced to B = 50mm, 
which is twice the fiber length. The length is Xmax = 800mm and the in- 
plane Eulerian mesh size is 2.5mm x 2.5mm. The Eulerian mesh size in 
thickness direction (z) ranges from 0.6 mm in the region of the final part 
geometry to 2mm at the top. The velocity in normal direction is set to zero 
at all boundaries of the domain. The domain of the initial stack is filled with 
40000 fiber bundles with a length of L = 25mm, which are discretized by ten 
truss elements each. Bundles interact with each other and with the mold 
surfaces through contact forces, but tangential lubrication is disabled due to 
its prohibitive restriction on the stable time step during explicit integration. 
The initial temperature distribution is precomputed assuming 10s contact 
with the lower mold between stack placement and the upper mold arriving 
at the stack. The pre-computed temperature and Eulerian mesh profile are 
depicted in Figure 8.22. 


Temperature in °C 
150.0 


100.0 


50.0 
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Figure 8.22: Initial stack at the flow front. The colors indicate the initial temperature distribution 
in all elements that are initially filled with material. Both mold halves are heated to a fixed value 
of 145°C. 


The SMC paste is a compressible, non-Newtonian and non-isothermal mate- 


rial with properties listed in Table A.8. Both mold halves are rigid body shells 
(element size 2.5mm x 2.5mm) with constraint degrees of freedom and a 
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prescribed temperature of 145°C. The upper mold’s vertical motion is con- 
trolled by a VUAMP subroutine that mimics the press controller, as described 
in Section 4.1.5. 


The Eulerian domain is subdivided at the sensor positions to evaluate the 
pressure as average of the nodal pressures which are filled with material at 
any given time. The compression force of the reference point linked to the 
upper rigid body mold is smoothed with a Butterworth filter, to eliminate 
spurious noise that occurs during the explicit time integration. 


8.3.3 Results and discussion 
8.3.3.1 75% mold coverage 


The recorded press data for 75% mold coverage is depicted in Figure 8.23 in 
light colors. The data lines represent the mean of five repetitions and the 
small standard deviation is depicted as light colored area behind the lines. 
All data is temporally aligned with the switch to force control, such that ft < 0 
refers to the displacement controlled profile and f > 0 refers to the force 
controlled profile. The mold is closed with a constant velocity of 1 mms! 
until it is controlled by the constant force resulting in a flatter slope and 
finally reaches its final thickness. For reference, the final thickness of the 
part is indicated with a dashed line and shows that the press displacements 
recordings are well calibrated. However, the final parts are approximately 
0.25 mm thinner than the expected thickness from volumetric considerations 
accounting for the measured compressibility. This is most likely caused by 
loss of material through the mold gap due to the high pressure experienced 
in this configuration. The pressure sensors show the expected behavior for 
a plug-flow with the left most sensor reaching the highest pressures and a 
decrease in pressure towards the flow front. However, the pressure of the left 
most sensor might reach its maximum capability of 200 bar and therefore 
might underestimate the maximum pressure in the mold. The resulting plates 
are homogeneous without any flow marks. 
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Figure 8.23: Press data at 75% mold coverage - experiment (light colors) and DBS model (dark 
colors). 


The DBS model is used to simulate the compression molding process in the 
press rheometer with 75% mold coverage and the result is shown in Figure 
8.23 in dark colors. As described previously, the numerical model predicts 
a thicker part than produced in the experiments, which is likely caused by 
leakage in the experiments. The displacement slope reduces correctly as 
soon as the virtual press controller reaches the force controlled phase. The 
time until the switch agrees with the experiments and the compression force 
is controlled reliably at the prescribed level. The spread between pressure 
sensors is similar to the experimental results and the general pressure levels 
are predicted correctly. However, the exact values are not met because of the 
earlier completion of the filling process in the simulation due to leakage, the 
maximum sensor capability of 200 bar and idealizations in the model. 
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8.3.3.2 50% mold coverage 


The experimental results with 50% initial mold coverage are depicted in 
Figure 8.24 in light colors. The compression force features a hump at which 
the rate of compression force increase slows down for a moment before 
reaching the maximum force. This phenomenon can be likely attributed to a 
transition from a sticking friction to a hydrodynamic lubricated friction. The 
final thickness of the parts agrees with the measured mold deformation and 
is in line with the expected thickness from mass conservation considering 
compressibility. 
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Figure 8.24: Press data at 50% mold coverage - experiment (light colors) and DBS model (dark 
colors). 


The measured pressures exhibit a larger variance between individual plates, 
which is indicated by the colored background areas representing the differ- 
ence between the upper and lower quartile of results. Contrary to the results 
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of 75% mold coverage, the left most pressure sensor P; does not record the 
highest pressures and the molded plates show evenly sized zones with dis- 
turbed fiber bundles exhibiting higher curvature at both ends of the molds 
(see Figure 8.25). The separation lines between these zones are displayed as 
gray lines in the mold coverage sketch in Figure 8.24. 


Figure 8.25: Photos of exemplary plates molded with 50% initial mold coverage. Both ends 
feature a darker zone with highly curved compressed fibers. 


The results of the DBS simulation with 50% mold coverage are displayed in 
Figure 8.24 in dark colors. The time to reach the compression force switch 
over point and the total compression time to fill the mold agree with the 
experimental data. The initial compression force increase and plateau are 
not met by the model even after accounting for an initial sticking behavior in 
the friction model. The used wall stress of 150 kPa should be sufficiently large 
to cause a steep compression force increase in the beginning according to 
the previous parametric study with one-dimensional model. An evaluation 
of the tangential mold stress reveals that this tangential stress is not properly 
enforced by Simulia Abaqus/Explicit at the reconstructed surface between 
molds and SMC, because it is only enforced if a surface overclosure is present. 
However, the contact initially oscillates between a closed and open state in 
many regions and thus allows a stack motion, whenever the contact is open. 
Hence, the initial compression force and the initial pressures are underes- 
timated. The DBS does not explicitly account for individual sheets and is 
therefore not capable to replicate the shear induced stack deformation that 
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causes the pressure sensor P; to report pressures below the values obtained 
by Po. 


8.3.3.3 25% mold coverage 


The experimental results for 25% initial mold coverage are depicted in Figure 
8.26 in light colors. This data was recorded on a smaller press, which deforms 
noticeably at the given load. Hence, the plot distinguishes the gap computed 
from press displacement and a true gap, which is corrected by the stiffness of 
mold and press. The true gap agrees with the measured final part thickness. 
The compression force increases in a relatively flat slope until the mold is 
almost filled and then raises quickly to the maximum compression force. 
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Figure 8.26: Press data at 25% mold coverage - experiment (light colors) and DBS model (dark 
colors). 
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The measured pressure does not adhere to the expected behavior of an ideal 
plug-flow (compare Figure 4.5), as the pressure does not monotonically de- 
crease from the wall towards the mold front. Instead, the leftmost pressure 
sensor P, records small pressures and at times the sensor Ps records the 
highest pressure. This is corroborated by the formation of large zones with 
disturbed fiber bundles, which are indicated by gray lines in Figure 8.26. De- 
spite the asymmetrical initial stack placement, these zones are approximately 
even in size at both ends of the mold. 


The formation of the zones is likely caused by a shear deformation of the stack. 
To analyze this effect, a sheet at the center of the nine-layered stack is colored 
by black paint to track its deformation. The black sheet is clearly visible after 
molding through the other sheets (see Figure 8.27) and its shape is similar to 
the central zone observed without coloring. This suggests that the central 
sheets are stretched less than predicted by the ideal plug-flow model and 
the upper and lower sheets fill the room at both ends, which originates from 
the reduced stretch of the central sheets. This mechanism is in line with the 
observations by Hohberg [158], but it is a stark contrast to the assumptions of 
the classical plug-flow models. Temperature gradients must affect more than 
just a very thin lubrication layer at the mold surface. 


Figure 8.27: Photos of exemplary plates molded with 25% initial mold coverage. The left plate is 
molded from nine sheets, were the central sheet was colored black. 


The DBS result with 25% initial mold coverage is displayed in Figure 8.26 in 
dark colors. The experimental results of this case were obtained on a smaller 
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press, which is compliant under the given loads and results in 0.2 mm appar- 
ently measured displacement at 4400 kN, even if there is none. Hence, the 
experimental gap determined from press displacement is corrected for the 
compliance by an empty reference press stroke. To model the compression 
process accurately, the virtual press controller mimics the press compliance. 
The resulting compression time and compression force agree with the ex- 
perimental data, but similar to the case with 5096 initial mold coverage, the 
initial increase and pressure drop of the left-most sensors due to the stack 
deformation are not captured. 


The orientation evolution in the plate is shown in Figure 8.28. The bundle ori- 
entations are mapped to an evaluation cell encompassing the entire domain 
using the procedure described in Section 5.1.6.3. The resulting orientation 
tensor represents the global orientation state in the cavity and its compo- 
nents are plotted as solid lines in Figure 8.28. The orientation has an initial 
preferred direction in the flow direction, because the removal of truss ele- 
ments at boundaries during the cutting procedure affects this direction less 
than the transverse direction with its longer boundary. Fiber bundles orient 
themselves in flow direction during the process, as expected. The solution 
of the one-dimensional reference model based on Jeffery's equation and the 
IBOF closure approximation is plotted as dotted lined for comparison. This 
solution assumes an ideal transversal isotropic initial fiber orientation state 
and is therefore shifted in comparison to the DBS result. However, except 
for the shift, this simple model agrees well with the mesoscale model. This 
suggests that the orientation evolution in a sufficiently large planar SMC 
region can be described equivalently with this macroscopic model and no 
further interaction parameters. 


8.3.4 Concluding remarks 
Press rheometer trials with distributed in-mold pressure sensors allow for an 
in-depth analysis of the compression molding process that allows to distin- 


guish contributions from compressibility, friction and the elongation of SMC 
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Figure 8.28: Orientation evolution for the entire cavity computed by DBS (solid) and with the 
macroscopic model using Jeffery's model (dotted). 


itself. The experimental results show that the SMC under investigation does 
not follow the ideal plug flow assumption for thick initial stacks. Contrary, 
the introduction of a colored sheet in the center of the stack proves that the 
central sheet is not stretched to the entire mold length and the outer layers 
fill the remaining region with a disturbed bundle structure. The bundle archi- 
tecture is not just a question of whether it is the region of the initial stack, but 
the result of complex flow mechanisms that can transport entire sheets to the 
center of the mold. 


The novel mesoscale simulation method can predict the compression time 
and occurring compression forces accurately with viscous behavior based 
only on a paste viscosity characterization at coupon scale. This coupon scale 
characterization is much simpler than press rheometry, but has the disadvan- 
tage that the pure paste might not be available for commercial SMC prepregs. 
The contributions from friction and anisotropic viscosity of the bundle sus- 
pension are proportionate. The non-isothermal simulation model accounts 
for compressibility of SMC, non-Newtonian viscosity, bundle-matrix interac- 
tions, bundle-bundle interactions (excluding lubrication), press compliance 
and press control. The predicted pressure differences between sensor po- 
sitions and overall pressure levels are predicted correctly, even if the exact 
pressures of the experiment are not met exactly. Reasons for the difference are 
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on the one hand experimental deficiencies, such as the maximum pressure 
capability of the sensors, leakage through the mold gap or the point-wise 
nature of the measurement, as well as on the other hand limitation in the 
simulation model, such as no explicit resolution of individual sheets and 
imperfect contact enforcement of the Simulia Abaqus/Explicit solver. Even 
though the initial sticking behavior cannot be accurately enforce in the con- 
tact, it is concluded that such a sticking friction is needed to model the entire 
compression process including squish and flow phase. Modeling the pres- 
sure drop of the first sensor and the motion of the inner sheet likely requires 
different fluid phases that represent individual SMC sheets. However, this 
would require a much finer Eulerian mesh and a solution for merging these 
phases with progressing simulation time, as heating and deformation remove 
any distinguishability of the sheets. 


8.4 Complex small part and co-molding 


A successor to the double curved dome in Section 8.1 has been designed 
in cooperation with Lucas Bretz and manufactured at KIT wbk, Karlsruhe, 
Germany. The part is used to validate the combined application of all DBS 
model features (non-isothermal, non-Newtonian, compressible, co-molding 
with patches) in a single part that contains beads and ribs. 


8.4.1 Experimental setup 


The part and initial stack position are shown in Figure 8.29. Two beads with 
different flank angles are positioned in the part and a rib is placed centrally 
in each bead. Both ribs are 2 mm wide at the top and have a draft angle of 3°. 
The rib height varies from 15 mm at one end to 0 mm at the other end. A lever 
is integrated into the part for easier recovery from the mold. The SMC stack 
has the initial dimensions 100 mm x 75 mm x 7 mm and is placed centrally 
into the mold, which is heated to 145 °C. Optionally, a single ply patch with 
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dimensions 80 mm x 7 mm x 0.45 mm can be placed in the narrow bead and 
a patch with dimensions 80 mm x 15 mm x 0.45 mm in the wide bead. The 
press profile is prescribed as a table with a maximum compression force of 
130 kN and parts are molded on a Lauffer hydraulic press in collaboration 
with Lucas Bretz at KIT wbk, Karlsruhe, Germany. 


120mm 


(a) Top view (b) Perspective view 
Figure 8.29: Complex small part geometry comprising two beads with different flank angles and 
centrally placed ribs with decreasing height. The molded part has outer dimensions 120 mm x 
94 mm. Patches are optionally placed inside the beads. 


The molded samples were analyzed by volumetric imaging using the same 
setup as described in Section 8.1, but with the resolutions 31 um and 50 um 
for the high resolution and low resolution scans, respectively. The uCT scans 
and volumetric image reconstruction were performed by Ludwig Schóttl at 
KIT IAM-WK, Karlsruhe, Germany. 


8.4.2 Numerical setup 


The process is modeled with the Direct Bundle Simulation technique ex- 
plained in Chapter 5. The upper and lower mold surfaces are extracted from 
the CAD model of the part, placed at an initial distance of 22.4 mm and are 
meshed with rigid body elements of 1 mm average edge length. The lower 
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mold remains fixed, while the upper mold is controlled by a virtual press con- 
troller (see 4.1.5) following the press profile listed in Table A.9. The Eulerian 
domain is expanded upward and downward beyond the mold surfaces and 
meshed with 0.8 mm x 0.8 mm x 0.8 mm EC3D8RT. The in-plane boundaries 
are implemented by constraining the velocity in normal direction to zero. The 
regions of initial stacks are filled with material in the beginning and 18000 
bundles of 25 mm are generated in the stack to represent the desired fiber 
volume fraction. The paste material is compressible, non-isothermal and 
non-Newtonian with parameterization from Section 7.1. A hydrodynamic 
friction model with parameters from Section 8.3 is employed at the mold 
walls. All simulation parameters are summarized in Table A.9. 


The basic configuration does not include unidirectional reinforcement 
patches. But patches can be added via 3.5 mm x 4mm and 3.75 mm x 4mm 
SART shell elements following the model described in Chapter 6 and with the 
parameterization in Section 7.2. The interaction between the patch surfaces 
and the SMC paste is modeled via no-slip conditions, as both utilize the same 
sticky matrix material. The interaction between the molds and patches is 
modeled with a Coulomb contact with a friction coefficient u = 0.3 according 
to the measurements by Bücheler [169]. 


8.4.3 Results and discussion 
8.4.3.1 Flow front progression 


The deformation of the SMC stack is shown in Figure 8.30 for the simulation 
and short-shots at corresponding thicknesses. First, the stack is formed to 
the shape of the mold without a flow process. As soon as the mold gap 
decreases below the thickness of the formed stack, a flow process starts and 
fills the ribs. It is noticeable that fiber bundles are pressed against the stamps, 
because both ends of fibers oriented in x direction are pulled in different 
directions, while the stamps press them down. Both ribs are filled equally 
from the bottom and the side connected to the bead wall. Hence, the last 
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filled region is not equivalent to the maximum rib height, but the volume 
farthest away from any entry to a rib. This behavior is coherent between 
simulation and experimental short-shots. As the mold is completely filled, 
a fraction of the simulated paste penetrates the rigid body walls and floats 
through the domain as droplets. 


Figure 8.30: Snapshots of the compression process with two beads and ribs. The left column 
shows the simulation results and the right column shows short shots of the molded parts. 
The complex geometry leads to leakage of the paste through the rigid mold surfaces due to 
inaccurate contact handling between the Eulerian surface and Lagrangian surface in Simulia 
Abaqus/Explicit. 
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8.4.3.2 Fiber-matrix separation 


Fiber-matrix separation can be observed with bare eyes in the molded parts 
due to the opaque resin. A exemplary part is scanned via uCT and a section 
across the rib in the wider bead is shown in Figure 8.31. Figure 8.31a shows a 
volumetric image reconstruction of the fiber architecture in the wider bead. 
The fiber bundles displayed here stay intact even after entering a rib. The 
slice of a lower resolution scan shown in Figure 8.31c reveals that a resin- 
rich region forms in the region that is filled last according to the simulated 
flow front progression. A higher resolution slice of this region is shown in 
Figure 8.31b. The mesoscale simulation approach predicts a similar pocket 
of resin in the same region, even though the actual detailed realization of 
the fiber architecture depends on the randomness of the initial charge in the 
simulation as well as the experiment. 


(a) 3D uCT view (b) High resolution uCT slice 


(d) Direct Bundle Simulation 
Figure 8.31: Fiber-matrix separation in the ribs. The uCT images were obtained by Ludwig 
Schóttl at KIT IAM-WK, Karlsruhe, Germany. 
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8.4.3.3 Patch damage 


The quality of parts molded with single patch layers varies significantly be- 
tween parts with intact patches and parts with defects such as fracture, dis- 
placement and deformation of patches. This is shown in Figure 8.32 and 
the variance is in line with the large scatter observed in the characterization 
experiments in Section 7.2. 


Figure 8.32: Patch defects observed in molded parts. 
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Figure 8.33: Simulation results of the section points adjacent to the mold surface. 


The result of a simulation with patches and fiber bundles is shown in Figure 
8.33. The flow front progress and resulting fiber architecture is unaffected 
by the presence of patches, except for a slightly faster filling process due to 
the volume occupied by patches in the cavity. The maximum temperature 
of patches adjacent to the mold surface is approximately 75 °C and therefore 
within the temperature range investigated during characterization. The patch 
model does predict damage and deformation in the patches, but the defor- 
mation mode does not strictly match the experimental results. Some paste 
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material reaches between patch and mold surface, because leakage at convex 
mold surfaces allows material to get to the other side of the patch. 


8.4.4 Concluding remarks 


The application example demonstrates the ability of the developed models 
to describe non isothermal compression molding and co-molding at com- 
ponent scale. The mold filling process and formation of resin-rich areas is 
predicted qualitatively correct and the uCT images illustrate that fiber bun- 
dles stay intact even after entering a rib. The partial penetration of Eulerian 
material through mold surfaces is a well known problem in compression 
molding simulations with the Coupled Lagrangian Eulerian framework of 
Simulia Abaqus/Explicit [158, 210]. It can occur at convex curved mold sur- 
faces due to the employed surface reconstruction method [194] of the internal 
solver. Large pressures and low viscosity foster the formation of droplets that 
penetrate the wall. Even though a finer mesh and smaller time step can 
reduce the effect, it always remains a disadvantage of the employed sur- 
face reconstruction technique. The integration of patches in the mesoscale 
compression molding process simulation is possible and the occurrence of 
damage in patches is predicted correctly. However, the accurate prediction 
of patch defects is difficult due to the inhomogeneous patch quality, which 
results in different experimental outcomes and uncertain model parameters. 


8.5 Demonstrator compression molding 


The demonstrator component of IRTG GRK 2078 serves as a use case that 
validates the applicability of the developed mesoscale DBS at larger scale 
components. The component is designed collaboratively within IRTG GRK 
2078 to incorporate three beads that reinforce the part under bending load. 
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8.5.1 Experimental setup 


The mold consists of a base mold with dimensions 800 mm x 250 mm and an 
insert geometry with dimensions 230 mm x 230 mm. This insert geometry is 
visualized in Figure 8.34a. The upper insert can be replaced with a ribbed 
structure that has varying rib features in each bead (see Figure 8.34b). The 
outer beads have six crossing rib pairs each, with a thickness of 3 mm and 
2mm, respectively. The radius at the rib base varies between 0.5 mm and 
1.5mm in these outer beads. The central bead has a higher rib density with a 
total of eight crossing rib pairs, but constant radii of 1.5 mm and a constant 
thickness of 2 mm. 


(a) Insert without ribs and insert-centered stack 


(b) Insert with ribs and full mold coverage 
Figure 8.34: Molds and initial stack configurations of the IRTG demonstrator component. 


The components are manufactured from UPPH-GF by Sergej Ilinzeer at Fraun- 
hofer ICT, Pfinztal, Germany on a Dieffenbacher COMPRESS PLUS DCP-G 
3600/3200 AS and a Dieffenbacher DYL630/500 hydraulic press with parallel 
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cylinder control. Compression force and cylinder displacements are recorded 
during these manufacturing molding trials. Parts without ribs are molded 
from a 320 mm x 240 mm x 7 mm initial stack that is placed centrally over 
the insert in the mold and parts with ribs are manufactured with full mold 
coverage of two SMC sheets. 


The demonstrator component' size prohibits an efficient evaluation of the 
fiber orientation via uCT for the entire component. However, the opacity of 
the employed resin system enables fluoroscopy with visible light. Hence, a 
fluoroscopy workflow was developed in collaboration with Sven Revfi from 
KIT IPEK, Karlsruhe, Germany and the workflow is depicted in Figure 8.35. 


g Camera 
~ 


Specimen 


Overhead Merge and 
projector processing 


Orientation, 
evaluation 


Figure 8.35: Image acquisition and image processing workflow. The specimens are placed on an 
overhead projector and images of the shine through region are taken. The images are merged and 
the lighting is homogenized in an image processing software to obtain a gray valued image of the 
entire specimen with 12.5 pixels mm”. The gray valued image is then analyzed with orientation 
analysis software. (The procedure was developed and applied together with Sven Revfi from KIT 
IPEK, Karlsruhe, Germany and Ludwig Schóttl from KIT IAM-WK, Karlsruhe, Germany.) 


Specimens are placed on an overhead projector and partial images are ac- 
quired with a Canon EOS 70D DSLR due to size restrictions and to improve 
the image resolution. These partial images are digitally processed to correct 
inhomogeneous lighting, distortions and improve contrast. Subsequently 
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they are aligned and merged to a single gray value image of the entire compo- 
nent. The procedure is repeated for the top and bottom side of three parts in 
total. Each image is analyzed with Orientation] [208] utilizing a smoothing 
length of 1 pixel and a grid size of 5 pixels at which orientation vectors are 
determined. In addition, the orientation on image data is analyzed using 
structure tensors [211]. The fluoroscopy procedure is corroborated by uCT 
with 42 mm resolution and 56 um resolution. The uCT scans, volumetric im- 
age reconstruction and fiber orientation analysis by structure tensors were 
performed by Ludwig Schóttl at KIT IAM-WK, Karlsruhe, Germany. 


8.5.2 Numerical setup 


The production process is modeled with the Direct Bundle Simulation tech- 
nique explained in Chapter 5. The upper and lower mold surfaces are ex- 
tracted from the CAD model of the part, placed at an initial distance of 
25.2 mm and are meshed with rigid body elements of 2.5 mm average edge 
length in the planar region and 1.0 mm average edge length in the bead region. 
The lower mold remains fixed, while the upper mold is controlled by a virtual 
press controller (see 4.1.5) following the press profile listed in Table A.10. The 
Eulerian domain encompasses the entire region between the molds and is 
meshed with 2.5 mm x 2.5 mm x 1.0 mm EC3D8RT elements in the planar 
regions and 1.0mm x 1.0 mm x 1.0 mm EC3D8RT elements in the region of 
the beads. The normal velocity of the boundaries of the Eulerian domain is 
constrained to zero and the region of the initial stack is filled with material 
in the beginning. For the unribbed components, a total amount of 150000 
bundles with L = 25mm length with element size / = 2.5 mm are generated in 
the stack to represent the desired fiber volume fraction. As the ribbed part is 
molded from a full initial mold coverage and the plate region has only little 
effect on the rib filling, the simulated domain is reduced to the insert area for 
this case. This also reduces the number of simulated fiber bundles to 30000, 
which are resolved finer with / = 1 mm. The paste material is compressible, 
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non-isothermal and non-Newtonian following the parameterization of Sec- 
tion 7.1. A hydrodynamic friction model with parameters from Section 8.3 is 
employed at the mold walls. All simulation parameters are summarized in 
Table A.10. 


8.5.3 Results and discussion 
8.5.3.1 Flow front progression 


Snapshots of the filling process are depicted in Figure 8.36. Initially, the 
stamps push the stack into the beads and fiber bundles are pulled into the 
shape leading to a higher orientation perpendicular to the beads. The outer 
stack contour is affected by this pull-in and the stack experiences out-of-plane 
deformation. At 4.8s, the beads are fully filled and the flow front extends 
uniformly in both directions. However, leakage occurs at some convex mold 
surface areas similar to the previous example in Section 8.4. Due to the 
asymmetric placement on top of the insert, the flow front reaches one end 
of the mold first (7.2s) and the flow direction reverts in some regions to 
completely fill the mold. 


Flow marks with highly aligned fibers form at the ends of the beads during 
molding (see Figure 8.37a, where flow marks are highlighted with a light green 
color). For an insert-centered initial stack, these flow marks are symmetric 
(equal length at both ends of the beads). The mesoscale DBS is able to predict 
the formation of these flow marks as a result of bundles with restricted motion 
due to the beads. As shown in Figure 8.37c, fiber bundles at the end of the 
beads travel a shorter distance than those surrounding them, which leads to 
an alignment of fibers. An adjustment of the initial stack to the mold center 
results in longer flow marks that are restricted to one side of the part (see 
Figure 8.37b, where flow marks are highlighted with a light green color). The 
effect of this small shift by 6cm can be reproduced by DBS, where similar 
long flow marks are predicted (see Figure 8.37d). The fiber architecture shows 
artifacts in the region of the beads, which are caused by leakage of paste 
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t=2.4s 


t=4.8s 


Figure 8.36: Snapshots of the simulated mold filling process for the IRTG demonstrator com- 
ponent without ribs and an insert-centered stack. The left column shows the top view on the 
component and the right column shows a section through the center bead. The last snapshot 
shows some droplets leaving the cavity due to leakage. 


droplets through the mold surfaces. These artifacts can occur, if a droplet 
leaves the domain rapidly and affects the neighborhood of a fiber bundle, 
with which it interacts. 
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(a) Insert-centered stack (photo) (b) Mold-centered stack (photo) 


(c) Insert-centered stack (DBS) (d) Mold-centered stack (DBS) 
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Figure 8.37: Formation of flow marks in the demonstrator part without ribs. The displacement 
describes the distance that each bundle traveled relative to its initial position. The insert- 
centered stack results in symmetric flow marks at all ends of beads. The mold-centered stack 
results in long one-sided flow marks, which are observed in experiments and the simulation. 


The ribbed demonstrator component molded with full mold coverage shows 
severe fiber-matrix separation in the ribs with resin-rich regions at the cross- 
ing points of each rib pair. The separation can be observed visually, as de- 
picted in Figure 8.38a, and seems to occur independent of rib width and radii 
at the rib entrance. The ribs are completely filled with paste in the mesoscale 
simulation, but fiber bundles are also segregated in the ribs, as shown in Fig- 
ure 8.38b. The segregation is caused by entangled bundles as well as bundles 
that have a limited capability to enter multiple ribs due to the inextensibility. 
Additionally, the drag force on bundle ends that just reach into a rib is not 
sufficient to move the entire bundle into the rib and increases separation. 


The central bead of a exemplary ribbed part is scanned by uCT and depicted 
in Figure 8.39a. The contour surface separating fiber bundles and matrix 
agrees with the visually observed separation. The DBS also shows a clear lack 
of bundles in the crossing points of ribs, as depicted in Figure 8.39b. 
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(a) Photo (b) DBS 
Figure 8.38: Overview on fiber-matrix separation in the insert region of the demonstrator com- 
ponent. (a) The dark translucent shade at the crossing points of ribs can be clearly attributed to 
fiber-matrix separation. (b) The crossing points of ribs are lacking fiber bundles. The matrix is 
not shown for visualization purposes, but the simulation predicts a complete filling of the ribs 
with matrix material and hence a fiber-matrix separation similar to the experiments. 


(a) uCT scan 


key as 


(b) DBS 
Figure 8.39: Overview on the fiber architecture of the central bead obtained from uCT and DBS. 


A more detailed view on the rib pair highlighted by a circle in Figure 8.39a is 
displayed in Figure 8.40 and compared to the corresponding simulation result. 
The flow front progression towards the center of the bead and resulting fiber 
architecture is predicted correctly. This is remarkable, as a recent study on 
very similar structures with commercial macroscopic models falsely predicted 
a preferred filling of the wider central crossing points of ribs [212]. The 
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authors assume that this stems from a lack of anisotropy in their models, 
but the results presented here suggest that the mesoscale structure plays an 
important role in the fill pattern of such ribs. 


(a) uCT scan (b) DBS 
Figure 8.40: Zoom into the fiber architecture in a central rib pair. The section is highlighted by a 
circle in Figure 8.38 


8.5.3.2 Fiber orientation 


The Aj; component of the fiber orientation tensor is depicted in Figure 8.41 
for the analysis with Orientation] on fluoroscopy images, structure tensors 
on fluoroscopy images, structure tensors on UCT scanned sections and the 
DBS result. The image data is evaluated on an exemplary sample and the 
DBS result is mapped on a 2.5mm x 2.5mm x 2mm evaluation mesh with 
the method described in Section 5.1.6.3. 


Both results based on image data show a predominant orientation in x- 
direction with a major perpendicular orientation in the beads and at the 
ends of the mold. This perpendicular orientation is more prominent in the 
structure tensor based computation. It is likely that this structure tensor 
based computation is more accurate due to better image pre-processing. 
Additionally, the structure tensor based computation leaves out the region 
with the sample label to avoid misinterpretation of the orientation in this 
region. Compared to image data, the evaluation of regions analyzed via uCT 
show a higher orientation in terms of absolute values, but the same trend of 
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(c) uCT (structure tensor) (d) DBS (see Section 5.1.6.3) 
HEU | o0 —EEEEEN 
0.0 0.5 1.0 
An 


Figure 8.41: Orientation tensor component Aj; for an exemplary part analyzed with all de- 
scribed methods. The computation in (b) leaves out the region with the sample label to avoid 
misinterpretation of the orientation in this region. 


perpendicular orientation in the beads. As the uCT scan contains informa- 
tion through the thickness of the sample, this suggests that the fluoroscopy 
method underestimates strong alignment due to its thickness averaging prop- 
erty. The DBS features distinct highly oriented regions at the end of each 
bead that correspond to the flow marks described in the previous section. 
These are less pronounced in the experimental data, even though strong fiber 
alignment is visually observed in the molded parts. The mesoscale simulation 
does predict a strong perpendicular orientation in the beads and even be- 
tween beads. This is reasonable, as fibers are pulled into the beads such that 
the initial stack is stretched ~140% in y-direction at the central area. Pulling 
fibers into the beads affects also the region between beads due to the fiber 
length of 25 mm. However, the effect is less pronounced in the experimental 
data and only visible in the central section of the uCT scan. More uCT scans 
of the other beads would be necessary to investigate the fiber architecture of 
this region in more detail. 
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8.5.3.3 Compression force 


The compression force of the demonstrator part is plotted in Figure 8.42. The 
plot shows the mean force of five trials in dark gray as well as the area between 
the lower quartile and upper quartile in light gray. The compression force 
features a plateau just after the initial forming phase, in which the stamps 
deform the stack to the beads. The most likely reason for this plateau is the 
initial sticking friction at the mold surface. The DBS simulation features a 
similar plateau, but it is less pronounced. This is likely affected by the poor 
contact enforcement in Simulia Abaqus/Explicit between the Eulerian phase 
and complex Lagrangian molds. Accordingly, the final thickness of the part is 
also underestimated due to the leakage at convex mold features. 
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Figure 8.42: Compression force and mold displacement during the demonstrator compression 
molding process. Dots mark the snapshot positions depicted in Figure 8.36. 


8.5.4 Concluding remarks 
The mesoscale Direct Bundle Simulation approach is applied to a large-scale 


component of dimensions 800 mm x 250 mm with ribs and beads. The model 
is able to predict molding defects such as the formation of flow marks and 
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fiber-matrix separation in ribs. In contrast to recent reports on macroscopic 
models in similar structures [212], the rib filling process is described qual- 
itatively correct from the outer borders towards the center and not with a 
preferred flow in the wider center of crossing rib pairs. The comparison of 
fiber orientation with experimental results shows a correct prediction of the 
predominant fiber orientation in x-direction of the plate with perpendicular 
orientation in the beads and at the end of the mold. However, the coarse 
mesh and previously described penetration of droplets through the mold 
walls cause artifacts in the bundle architecture and an underestimation of 
the final part thickness. The mapping procedure of the orientation evalua- 
tion can be also used to transfer the process simulation results to a mesh for 
subsequent structural simulation. Contrary to regular macroscopic models 
applied in the state of research, the new mesoscale model allows to transfer 
fiber volume fraction and evaluate higher order fiber orientation tensors, 
if desired. A comparison to results obtained with the commercial software 
Autodesk Moldflow may be found in the thesis of Revfi [213]. 


210 


9 Summary 


After a summary of the state of research in SMC compression molding and 
stating the objective, this work proposes a macroscale reference model for 
SMC process simulation. The model is supplemented by utility functions 
such as a virtual press controller and a friction model, which are also used for 
the mesoscale simulations. The main contribution is the development of a 
mesoscale Direct Bundle Simulation (DBS) approach that models SMC paste 
as a combination of paste moving in an Eulerian frame and Lagrangian truss 
elements representing individual fiber bundles. The model includes bundle- 
bundle interactions as well as bundle-matrix interactions. A collection of 
pre- and post-processing tools for productive application at component scale 
completes the mesoscale model. Additionally, an anisotropic patch model 
with individual damage variables for stitching yarn and infiltrated UD carbon 
fibers is developed. The required parameters of the models are characterized 
and applied to several validation examples. 


9.1 Conclusions 


The prediction of fiber architecture via the novel mesoscale Direct Bundle 
Simulation is improved in comparison to macroscale models in confined 
regions and close to the mold walls. The method is able to predict the for- 
mation of knit lines and flow marks in the same positions as experimentally 
validated. All this can be applied for entire components and is to the au- 
thors best knowledge the first method resolving individual bundles enabling 
component simulations. 
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For large planar SMC regions, Jeffery’s model without interaction parameters 
seems to be a good choice for a macro model, as it agrees with the more 
detailed DBS in such large planar regions. This is remarkable, because it 
does not require additional interactions parameters and retarding principal 
rate parameters as they have been developed for injection molding. The 
macroscopic model features up to ten times faster computation times than 
the mesoscale model, but it is not able to predict fiber-matrix separation and 
bundle curvatures. 


The proposed mesoscale method allows the prediction of fiber-matrix sepa- 
ration and bundle curvature, because flexible bundles and paste are able to 
move at different speeds. The predicted locations of resin-rich regions agree 
with experimental results. 


Lubrication forces between individual fiber bundles are incorporated into 
the model, but are numerically challenging due to the resulting small time 
step during explicit integration. The introduction of lubrication affects the 
compression force during molding of planar SMC parts, while the orientation 
is hardly affected. In most cases, mold friction is a dominant contributor 
to compression force in comparison to lubrication and lubrication may be 
neglected in favor of faster computation times. 


The ideal plug-flow assumption has limited validity for the considered UPPH- 
GF SMC and sheets do slip relative to each other, which has been proven by a 
colored sheet. Even with mesoscale DBS and accounting for mold friction, 
press compliance, compressibility, non-Newtonian non-isothermal paste, it 
is not possible to describe this layer-induced mechanism with a single matrix 
phase. Likely, it would require separate matrix phases, which arises new 
challenges w.r.t mesh resolution and the ability to merge such phases. 


To model the squish and flow phase of SMC, sticking friction should be added 
to the hydrodynamic friction model. Therefore, a friction model is suggested 
that transitions from an initial sticking behavior to a hydrodynamic power- 
law friction model. However, the friction model relies on the implementation 
of the surface reconstruction from the Eulerian phase. The employed Simulia 
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Abaqus/ Explicit CEL method shows some weaknesses in this area that lead to 
an imperfect contact and leakage of paste material through the mold surfaces. 


Patches from the in-house production suffer from significant scatter due to 
a missing series production quality control. They can be integrated to the 
simulation and the occurrence of damage is predicted correctly, but the exact 
deformation mode cannot be predicted reliably due to the scatter in patch 
properties. 


9.2 Outlook and recommendations 


For macroscopic SMC orientation models, the mesoscale simulation results 
suggest that a focus on further development of refined fiber orientation 
parameters is not necessary for large planar SMC regions. Instead, it seems 
more promising to further develop models that account for confinement in 
narrow regions following for example the work of Perez et al. [25]. 


The mesoscale process model can enrich the previously developed CAE chain 
for SMC [6] with higher-order fiber orientation tensors, fiber volume frac- 
tion, bundle curvature and other information given by the more detailed 
prediction of fiber bundle architecture. Given the probabilistic nature of 
some molding mechanisms, the scatter should be assessed with a continuous 
probabilistic CAE chain that leverages the ability of the developed model to 
predict different realizations for different initial randomly generated stacks 
under identical processing conditions. 


The mesoscale model itself could be improved by extensions to model flat- 
tening of bundles and disintegration of bundles that can occur in rare cases. 
For example, a damage variable could describe the disintegration of the bun- 
dle character based on the subjected shear stress on bundles. Finally, the 
use of an improved implementation could address the leakage problem and 
generally improve the computation performance in Simulia Abaqus/ Explicit 
CEL. 
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A Appendix 


A.1 Sketches for geometric considerations 


Equation (5.20) was obtained by a projection, which is illustrated in Figure 
A.la. The normal direction is computed from the distance between the tip 
of pj and the projection of p; on the unit vector in direction of the velocity 
[[Av;]]. 


The curvature given in Equation (7.4) is the inverse radius of the bent section. 
The radius Rp comprises a lower section with length ro sin ay and an upper 
section with length (ro + ro cos(ap)) cota. 


dj 


(a) Normal direction (b) Bending 
Figure A.1: Sketches for the derivation of Equation (5.20) and Equation (7.4) 
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A.2 Simulation parameters 


Table A.1: Parameter variation of the reference model presented in Figure 4.6 in Section 4.3. 


Description Symbol Default value Changed value 
Initial length Xo 0.29m 
Width Bo 0.25m 
Density po 1900 kg m? 
Initial orientation Axx,0 0.5 

Ayyo 0.5 

Axy,0 0.0 
Fiber volume fraction f 0.0 
Shear viscosity 7 1000 Pas 2000 Pas 
Bulk modulus Ko 190 000 Pa 1900 Pa 
Short-range interaction kp 0.0 
Shear factor Ns 0.0 
Friction A 0Nsm^? 50Nsm^? 
Force Fmax 1000 
Profile ho, ho 9mm, -1 mm st 

hy, hy 0mm, -1 mms 
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A.2 Simulation parameters 


Table A.2: Simulation parameters for anisotropic flow verification in Section 5.2.2. 


Parameter Symbol Isotropic Macro DBS 

Truss element size 1 - - 2.5mm 

Fluid element size - 2.5 mm x 2.5 mm x 1.0 mm 

Stack dimensions - 25 mm x 25 mm x 15 mm 

Density Po 1900 kg m? 

Initial orientation Axx,0 OS 0.0 0.0 
Ayyo 0.5 1.0 1.0 
Axyo 0.0 0.0 0.0 

Fiber volume fraction f 0.25%-2.5% 

Matrix viscosity n 0.1Pas 

Aspect ratio Tp 50 

Nominal bundlelength L 10mm 

Bulk modulus Ko 1500 Pa 

Short-range interaction kp 0 0 - 

Shear factor Ns 0 0 - 

Friction A 0Nsm^? 

Profile ho, ho 20mm, -1 mm sl 
hy, hy 0mm, -1 mms! 


Table A.3: Simulation parameters for linear compression example with 0% fiber volume fraction. 


Parameter Symbol 1D model 3D model DBS 
Initial length Xo 50mm 
Width Bo 50mm 
Fluid element size - 1.25mm 1.0mm x 1.0 mm x 1.0mm 
Density Po 1900 kg m? 
Mass scaling Km 1000 
Matrix viscosity 1 100kPas 
Bulk modulus Ko 190 MPa 
Friction A 0Nsm^? 
Profile ho, ho 5mm, -1.667 mms! 
hy, hy 0mm, -1.667 mm s 
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Table A.4: Simulation parameters for linear compression example with 1% fiber volume fraction. 


Parameter Symbol 1D model 3D model DBS 
Initial length Xo 50mm 
Width Bo 50mm 
Truss element size 1 - - 2.5mm 
Fluid element size - 1.25mm 1.0 mmx 1.0 mm x 1.0 mm 
Density Po 1900 kg m? 
Mass scaling Km 10000 
Initial orientation Axx,0 0.5 0.5 0:5 
Ayyo 0.5 0.5 0.5 
Axy,0 0.0 0.0 0.0 
Fiber volume fraction f 1.0% 
Aspect ratio Tp TAL) 122) - 
Nominal bundle length L - - 25mm 
Equivalent diameter d - - 0.2mm 
Matrix viscosity n 100kPas 
Bulk modulus Ko - 190 MPa 190 MPa 
Bundle elastic modulus E - - 72 GPa 
Friction A 0Nsm^? 
Profile ho, ho 5mm, -1.667 mm st 
hy, hy 0mm, -1.667 mm s7} 
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A.2 Simulation parameters 


Table A.5: Simulation parameters for linear compression example with 25% fiber volume fraction. 


Parameter Symbol 1D model 3D model DBS 
Initial length Xo 50mm 
Width Bo 50mm 
Truss element size 1 - - 2.5mm 
Fluid element size - 1.25mm 1.0mmx1.0mmx1.0mm 
Density Po 1900 kg m? 
Mass scaling Km 10000 
Initial orientation Axx,0 0.5 0.5 0.5 
Ayy,0 0.5 0.5 0.5 
Axy,0 0.0 0.0 0.0 
Fiber volume fraction f 25.096 
Aspect ratio Tp 172,8) 72.8) - 
Nominal bundle length L - - 25mm 
Equivalent diameter d - - 0.2mm 
Matrix viscosity n 100kPas 
Bulk modulus Ko - 190MPa 190 MPa 
Bundle elastic modulus E - - 72 GPa 
Friction A 0Nsm^? 
Profile ho, ho 5mm, -1.667 mm s 
hy, hy 0mm, -1.667 mms! 


219 


A Appendix 


Table A.6: Simulation parameters for simple double curved dome. 


Parameter Symbol Macro DBS 
Truss element size il - 2.5 mm 
Fluid element size - 1.0 mmx 1.0 mm x 1.0 mm 
Density po 1900 kg m? 
Initial orientation Axx,0 0.5 0.5 
Ayy,o 0.5 0.5 
Axyo 0.0 0.0 
Fiber volume fraction f 23% 
Matrix viscosity n 1x 10° Pas 
Nominal bundle length L 25 mm 
Bulk modulus Ko 1.9 GPa 
Short-range interaction kp 0 - 
Shear factor Ns 0 - 
Bundle elastic modulus E - 72 GPa 
Bundle density Pb - 2600 kg m^? 
Bundle radius R - 0.1 mm 
Friction A 1 x 10 Ns m` 
0.6 
vo 0.001 ms! 
Mass scaling factor Km 1x 106 
Time step At 3x10*s 
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A.2 Simulation parameters 


Table A.7: Simulation parameters for the honeycomb part. 


Parameter 


Symbol Value 


Truss element size 
Fluid element size 
Density 

Initial orientation 


Fiber volume fraction 
Matrix viscosity 
Nominal bundle length 
Bulk modulus 

Bundle elastic modulus 
Bundle density 

Bundle radius 

Friction 


Mass scaling factor 
Max. compression force 
Initial gap 

Press profile (B) 


I 2.5mm 


- 1.0mm x 1.0 mm x 1.0 mm 


po 1480 kg m? 

Axx,0 0.5 

nmm 0.5 

Axyo 0.0 

f : 36.6% 

n 25000 Pas 

JL, 50mm 

Ko 1.48 GPa 

iE 73 GPa 

Pb 2600 kg m^? 

R 0.237 mm 

A 5.59 x 10 N s m~’ 

m 0.66 

vo 0.001 ms”! 

Km 1x 10% 

Finax 1500 kN 

ho 12mm 

hy, hy 14mm, -10mm s7} 

h2, h2 10 mm, -2.5 mm sl 
1 


ha, hg 0mm, -10 mm s” 
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Table A.8: Simulation parameters for the press rheometer. 


Parameter Symbol Value 
Truss element size l 2.5mm 
Fluid element size - 2.5mm x 2.5 mm x 0.6 mm 
Density po 1480 kg m3 
Initial orientation Axx,0 0.5 
Ayyo 0.5 
Axy,0 0.0 
Fiber volume fraction f 23% 
Matrix viscosity Di 72000 Pas 
yo 0.1 
n 0.385 
Jp 40.73°C 
Al 7.94 
A2 105.69 °C 


Thermal conductivity « 


Specific heat capacity cp 
Nominal bundle length L 


Equation of state 


Bundle elastic modulus E 


Bundle density Pb 
Bundle radius R 
Friction A 
vo 
TQ 
Ut 
Gap conductance kr 
Mass scaling factor Km 


Max. compression force Fmax 


Press profile 
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0.163Wm kl 
1530]kg ! K! 
25mm 
see Table 7.4 
73 GPa 
2600 kg m^? 
0.1 mm 
3.0 x 10 N sm? 
0.6 
0.001 ms! 
150kPa 
10 mms” 
403W m”? K! 
1x 10% 
4400 kN 
9mm,-1mms- 


1 


1 


0mm, -1mms' 


A.2 Simulation parameters 


Table A.9: Simulation parameters for the complex small part. 


Parameter Symbol Value 
Truss element size I 1mm 
Fluid element size - 1mmximmximm 
Density Po 1480 kg m? 
Initial orientation Axx,0 0.5 

Ayyo 0.5 

Axyo 0.0 
Initial temperature To 25°C 
Fiber volume fraction f 23% 
Matrix viscosity Dı 72000 Pas 

yo 0.1 

n 0.385 

TS 40.73°C 

Ay 7.94 

A» 105.69 °C 
Thermal conductivity « 0.163 Wm! K~! 
Specific heat capacity ^ cp 1530] kg! K-! 
Nominal bundle length L 25mm 
Equation of state p(Ae) see Table 7.4 
Bundle elasticmodulus E 73GPa 
Bundle density Pb 2600 kg m^? 
Bundle radius R 0.1 mm 
Friction A 3.0 x 108 N s m~’ 

m 0.6 

vo 0.001 ms! 

TQ 150kPa 

Ut 1l0mms! 
Gap conductance kr 403WmK! 
Mass scaling factor Km 1x 10° 
Max. compression force Fmax 130kN 
Initial gap ho 22.4mm 
Press profile hy, hy 25mm, -35 mm sl 

ho, ho 10mm, -25mms 

ha, h3 5mm, -15mm sl 

ha, ha 2.5mm, -3mms'! 

hs, hs 0mm, -3mm sl 
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Table A.10: Simulation parameters for the IRTG demonstrator part. 


Parameter Symbol Without ribs With ribs 
Truss element size Į 2.5mm 1.0mm 
Fluid element size - 2.5mm x 2.5 mm x 1 mm 1mmximmximm 
Density po 1480 kg m? 
Initial orientation Axx,0 0.5 
Ayy,o 0.5 
Axyo 0.0 
Initial temperature To 308€ 25°C 
Fiber volume fraction f 23% 18% 
Matrix viscosity Di 72000 Pas 
Yo 0.1 
n 0.385 
D 40.73°C 
ET 7.94 
Ao 105.69*C 
Thermal conductivity — x 0.163 Wm! K~! 
Specific heat capacity Cp 1530 Jkg! K! 
Nominal bundle length L 25mm 
Equation of state p(Ae) see Table 7.4 
Bundle elasticmodulus E 73GPa 
Bundle density Db 2600 kg m^? 
Bundle radius R 0.1 mm 
Friction A 3.0 x 10 N s m`’ 
0.6 
vo 0.001 m s7} 
To 150 kPa 
Vt 10mms'! 
Gap conductance kr 403Wm~ kK! 
Mass scaling factor Km 1x10° 
Max. compression force Fmax 1400 kN 
Initial gap ho 25.2mm 22.2mm 
Press profile hy, hy 50mm, -30 mm sl 23mm, -15mm sl 
ho, ho 10mm, -5mm s7} 5mm,-1mms'! 
h3, h3 5mm, -1mms 0mm, -1 mms! 
hs, hs 0mm, -1 mms”! 
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